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A longstanding limitation of first-principles calculations of substitutional alloy phase diagrams is 
the difficulty to account for lattice vibrations. A survey of the theoretical and experimental litera- 
ture seeking to quantify the impact of lattice vibrations on phase stability indicates that this effect 
can be substantial. Typical vibrational entropy differences between phases are of the order of 0.1 
to 0.2fcs/atom, which is comparable to the typical values of configurational entropy differences in 
binary alloys (at most 0.693fcs /atom). This paper describes the basic formalism underlying ab initio 
phase diagram calculations, along with the generalization required to account for lattice vibrations. 
We overview the various techniques allowing the theoretical calculation and the experimental de- 
termination of phonon dispersion curves and related thermodynamic quantities, such as vibrational 
entropy or free energy. A clear picture of the origin of vibrational entropy differences between phases 
in an alloy system is presented that goes beyond the traditional bond counting and volume change 
arguments. Vibrational entropy change can be attributed to the changes in chemical bond stiffness 
associated with the changes in bond length that take place during a phase transformation. This 
so-called "bond stiffness vs. bond length" interpretation both summarizes the key phenomenon 
driving vibrational entropy changes and provides a practical tool to model them. 
Submitted to Reviews of Modern Physics. 



I. INTRODUCTION 



The field of first-principles alloy theory has made substantial progress over the last two decades. It is now possible 
to predict relatively complex solid-state phase diagrams starting from the basic principles of quantum mechanics 
and statistical mechanics. Since no experimental inpu t is required , thes e ab-initio calculations h ave been useful for 
clarifying the phase diagram of several new materials (Gcder et al., 1990; ran dcr Vcn ct al, 199S). Several excellent 



reviews on the topic exist ( Ducastcllc , 1991 ; dc Fontaine , 1994 ; Zungcr , 1994]) . The accuracy of calculated phase 



diagrams is currently limited by two factors. First, one needs, as a starting point, the energy of the alloy in various 
atomic configurations and hence, one is limited by the accuracy of the quantum-mechanical calculations used to 
obtain these energies. Typically, methods based on Density Functional Theory (DFT), such as the Local Density 
Approximation (LDA) or the Generalized Gradient Approximation (GGA) are used. 

A second shortcoming arises from the fact that, in order to reduce computational requirements, the sampling of 
the partition function to obtain the free energy is only done over a limited number of degrees of freedom. Typically, 
these include substitutional interchanges of atoms but no atomic vibrations. Attempts to either assess the validity 
of this approximation or to devise computationally efficient ways to account for lattice vibrations are currently the 
focus of intense research. This interest is fueled by the observation that phase diagrams obtained from first principles 
often incorrectly predict transition temperatures. It is hoped that lattice vibrations could account for some of the 
remaining discrepancies between theoretical calculations and experimental measurements. 

Three main questions are addressed in this paper. 

1. Do lattice vibrations have a sufficiently important impact on phase stability that their thermodynamic effects 
need to be included in phase diagram calculations? 

2. What are the fundamental mechanisms that explain the relationship between the structure of a phase and its 
vibrational properties? 

3. How can the effect of lattice vibrations be modeled at a reasonable computational cost? 

This paper is organized as follows. First, Section |l| presents the basic formalism that allows the calculation of 
phase diagrams, along with the generalization needed to account for lattice vibrations. A review of the theoretical 
and experimental literature seeking to quantify the impact of lattice vibrations on phase stability is then presented 
in Section [II. The main mechanisms through which lattice vibrations influences phase stability are described in 
Section IV . Section [v] describes the methods used to calculate vibrational properties while Section VI presents the 
experimental techniques allowing their measurement. Finally, Section VII discusses the strengths and weaknesses of 
a variety of models of lattice vibrations. 
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II. THE FORMALISM OF ALLOY THEORY 



Phase stability at constant temperature is determined by the free energ yQ F. The free energy can be expressed as 
a sum of a configurational contribution F con fi g and vibrational contributions -F V ib- The configurational contribution 
accounts for the fact that atoms can jump from one lattice site to another, while vibrational contribution accounts 
for the vibrations of each atom around its equilibrium position. The first part of this section presents the traditional 
formalism used in alloy theory to determine the configurational contribution. The second part introduces the basic 
quantities that determine whether lattice vibrations have a significant effect on phase stability. The third part 
describes how the traditional formalism can be adapted when lattice vibrations do need to be accounted for. 



A. The cluster expansion formalism 



One of the goals of alloy theory is to determine the relative stability of phases characterized by a distinct ordering of 
atomic species on a given periodic array of sites. This array of sites, called the parent lattice, can be any crystallographic 
lattice augmented by any motif. A convenient representation of an alloy system is the Ising model. In the common 
case of a binary alloy system, the Ising model consists of assigning a spin-like occupation variable o~i to each site i of 
the parent lattice, which takes the value —1 or +1 depending on the type of atom occupying the site. A particular 
arrangement of spins of the parent lattice is called a configuration and can be represented by a vector a containing 
the value of the occupation varia ble for each sit e in th e parent lattice. Although this framework can be extended to 
arbitrary multicomponent alloys (Sanchez et al. , 1984), we focus on the case of binary alloys, since all the studies we 
review consider binary alloys only. 

When all the fluctuations in energy are assumed to arise solely from configurational change, the Ising model is a 
natural way to represent an alloy. The thermodynamics of the system can then be summarized in a partition function 
of the form: 



Z = J2^P(-PE(a)) (1) 

a 

where (3 = 1/ (ksT), and E{a) is the energy when the alloy has configuration a. It would be computationally 
intractable to compute the energy of every configuration from first-principles. Fortunately, the configurationa l depen- 
dence of the energy can be parametrized in a compact form with the help of the so-called cluster expansion ( [Sanchez] 



ct al., 1984). The cluster expansion is a generalization of the well-known Ising Hamiltonian. The energy (per atom) 



is represented as a polynomial in the occupation variables: 

E(a) 



N 




= ^2 m a J a ( Y[ <?i ) (2) 



where a is a cluster (a set of sites i). The sum is taken over all clusters a that are not equivalent by a symmetry 
operation of the space group of the parent lattice, while the average is taken over all clusters a' that are equivalent 
to a by symmetry. The coefficients J a in this expansion embody the information regarding the energetics of the alloy 
and are called the effective cluster interaction (ECI). The multiplicities m a indicate the number of clusters that are 
equivalent by symmetry to a (divided by the number of lattice sites) ^] 

It can be shown that when all clusters a are considered in the sum, the cluster expansion is able to represent any 
function E (a) of configuration a by an appropriate selection of the values of J a . However, the real advantage of 
the cluster expansion is that, in practice, it is found to converge rapidly. A sufficient accuracy for phase diagram 
calculations can be achieved by keeping only clusters a that are relatively compact (e.g. short-range pairs or small 
triplets) . The unknown parameters of the cluster expansion (the ECI) can then determined by fitting them to the 
energy of relatively small number of configurations obtained, for instance, through first-principles computations. The 



1 Strictly speaking, at constant pressure, the Gibbs free energy G = F + PV should be used instead of the Helmoltz free energy 
F, but at atmospheric pressure, the PV term is negligible for an alloy. 

2 Both the number of clusters and the number of sites are infinite but their finite ratio can be obtained by ignoring all but one 
periodic repetitions of the clusters (or the atoms) by the a translational symmetry operation of the lattice. 
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cluster expansion thus presents an extremely concise and practical way to model the configurational dependence of 
an alloy's energy. 

How many ECI and structures are needed in practice? A typical well-converged cluster expansion of the energy 
of an alloy consists of about 20 to 30 ECI and necess i tates the calculation of the e nergy of around 40 t o 50 ordered 



structures (see, for instance, (van dcr Vcn ct al., 1998; Garbulksy and Ccdci, 1995| ; pzolins ct al. , 1998c)). A faithful 



modeling of the qualitative features of the phase diagram (correct stable phases and topology) typically requires far 
fewer EC I (as little as 1 pair interaction ) and c orrespondingly less structures, as illustrated by the numerous examples 
given in ( de Fontaine , 1994 ; Ducastelie] , 1991). In general multicomponent systems, the number of ECI and ordered 
struct ures r equired to achieve a given precision unfortunately grows rapidly with the number of species ((Sanchez 
et a|], |1984|) ). For instance, in ternaries, each pair interaction is characterized by 3 interaction parameters instead of 
only one in the binary case. For this reason, v e ry few first- principle calculatio ns of ternary phase diagrams have been 
attempted (see (McCormack and dc Fontaine, 1996| ) for a recent example, or (dc Fontaine, 1994) for a survey). 

Although the cluster expansion usually allows a very compact representation of the energetics of an alloy system, 
there are two situations where a standard cluster expansion is known to converge slowly. Systems where long-range 
elastic interactions are important due to a large atomic size mismatch between the alloyed species may require 
t hat clastic interactions be explicitly accounted for t hrough the use of a so-called reciprocal space cluster e xpansion 



( ^ungei , 1994 ; Ozolins et al. , 19986 ; Ozolins et al. , 1998c ). Another situation, as recently identified by ( Johnson 



et a|J, 200C ) , is when the electronic structure of the system exhibits a very strong configurational dependence due to 
symmetry-breaking effects. 

In the cases where a short-range cluster expansion does provide a sufficient accuracy, the process of calculating the 
phase diagram of an alloy system can be summarized as follows. First, the energy of the alloy in a relatively small 
number of configurations is calculated, for instance through first-principles computations. Second, the calculated 
energies are used to fit the unknown coefficients of the cluster expansion (the ECI J a ). Finally, with the help of this 
compact representation, the energy of a large number of configurations is sampled, in order to deter mine the phase 
boundaries . This latter step can be accomplished with eit her the Cluster Var iation Method (CVM) (Kikuchi, |1951 



Ducastclle , 1991), th e low-temperature expansion (LTE) ( Kohan et al. , 1998 ), or Monte-Carlo simulations (Binder 



and 



Heermann , 



1988). 



B. The effect of lattice vibrations 



The previous Section described the framework allowing the calculation of phase diagrams under the assumption 
that the thermodynamics of the alloy is determined solely by configurational excitations. Accounting for vibrational 
excitations introduces corrections to this simplified treatment. This section presents the basic quantities that enable 
an estimation of the magnitude of the effect of lattice vibration on alloy thermodynamics. To understand the effect 
of lattice vibrations on phase stability, it is instructive to decompose the configurational ( "config" ) and vibrational 
( "vib" ) parts of the free energy F a of a phase a into an energetic contribution E and an entropic contribution S : 

r — ^config — 1 D coiifig "+" ^vib — 1 D vib- W 

We take the convention that E" onfig is the energy of the alloy system when all atoms are frozen at their average position 
at a given temperature. In the approximation of harmonic lattice vibrations and in the limit of high temperature, 
the vibrational energy E" ih is simply determined by the equipartition theorem and is independent of the phase a 
considered. Hence as long as these approximations are appropriate, lattice vibrations are mainly expected to influence 
phase stability through their entropic contribution S^ ih . 

Intuitively, the vibrational entropy S" ib is a measure of the average stiffness of an alloy, as can be best illustrated 
by considering an simple system made of large number of identical harmonic oscillators. The softer the oscillators 
are, the larger their oscillation amplitude can be, for a fixed average energy per oscillator. Hence, the system samples 
a larger number of states and the entropy of the system increases. In summary, the softer the alloy, the larger the 
vibrational entropy. 

A phase with a large vibrational entropy is stabilized relative to other phases, since a larger vibrational entropy 
results in a lower free energy, as seen by Equation (|3|) . From a statistical mechanics point of view, this fact can be 
understood by observing that a phase that encloses more states in phase space is more likely to be visited, as the 
system undergoes microscopic transitions, and therefore exhibits an increased stability. 

The central role of vibrational entropy can be further appreciated by considering the effect of vibrations on a 
phase transition between two phases a and (3 which differ only by their average configuration {e.g. an order-disorder 
transition). If the vibrational entropy difference between the two phases is AS"^ 13 , the transition temperature 
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obtained with both configurational and vibrational contributions (^ c "nfig+vib) i s related to the transition temperature 
obtained with configurational effects only (T c "^) by 



T 



config+vib 



t: 



as: 



config 



^ D config, 



(4) 



where AS^^fg is the change in configurational entropy upon phase transformation (Garbulsky and Cedei, 1994) 
This result is exact in the limit of small vibrational effects, hig h temperatur e and h armonic vibrations. A correction 



to this result that accounts for anharmonicity can be found in (Ozolins ct al., 1998c). Equation M) indicates that the 



quantity determining the magnitude of the effect of lattice vibration on phase stability is the ratio of the vibrational 
entropy difference to the configurational entropy difference. For this reason, most investigations aimed at assessing 
the importance of lattice vibrations focus on estimating vibrational entropy differences between phases. Since the 
configurational entropy (per atom) S CO nfig for a binary alloy at concentration c is bracketed by 

< Config < -k B (cine + (1 - c)ln(l - c)) < fc s ln2 » 0.693/c B , (5) 

Equations ([|) and (||) provide us with a absolute scale to gauge the importance of vibrations. As we will see, typical 
vibrational entropy differences are of the order of 0.2fcs , indicating that corrections of the order of 30% to the transition 
temperature may not be uncommon. 

While it is clear that vibrational excitations introduce quantitative corrections to the simple picture of alloy ther- 
modynamics based on configurational excitations only, more profound effects of a qualitative nature are also possible. 
Vibrational effects may lead to deviations f rom the traditional bel i ef tha t , at high eno ugh temperature, all short-range 
order in a disordered material disappears (Garbulsky and Ceder, 1994; Miller, 1989] ). While a fully disordered state 
clearly maximizes configurational entropy S'config, it is not clear that the total entropy S con fig + S v ib is necessarily 
maximized in the state of maximum configurational disorder. The presence of short-range order may increase the 
total entropy, relative to a fully disordered alloy, through an increase of the vibrational entropy. Vibrational entropy 
somewhat challenges our intuition, which is largely derived from tacitly assuming that configurational disorder is all 
disorder. It is even conceivable that vibrational entropy could induce a transition from a disordered to an ordered 
phase with increasing temperature, if the vibrational entropy difference between the ordered and disordered phases is 
larger and opposite to the configurational entropy difference. While this phenomenon has, so far, not been observed 
in metallic alloys, presumably because of the large configurational entropy associated with disordering, it does occur 
in molecular systems, such as in diblock copolymer melts, where the configurational entropy (per monomer) is small 
([Russell ct all |1994|). 



C. Coarse graining of the partition function 



The cluster expansion formalism presented in Section II A appears to focus solely on configurational excitations. 
This section, shows that, in fact, non-configurational sources of energy fluctuations can naturally be taken into 
account within the cluster expansion framework through a process called "coarse graining" of the partition function^] 



(Cedei, 1991; Ceder, 1993). This procedure also clarifies the nature of the physical states that are represented by a 
configuration of the Ising model. 

All the thermodynamic information of a system is contained in its partition function: 



Z = J2^P [-ffi], 



(6) 



where Ei is the energy of the system in state i. In the case of a crystalline alloy system, the sum over all possible 
states of the system can be conveniently factored as following: 



Z = EEEE ex P \-PE{L, a, v, e)] 



where 



(7) 



L cr£L v£a edv 



This "coarse graining" process is, of course, unrelated to any processing aimed at increasing the grain size in a polycrystallmc 
material. 
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• L is a so-called parent lattice: it is a set of sites where atoms can sit. In principle, the sum would be taken over 
any Bravais lattice augmented by any motif. 



• a is a configuration on the parent lattice: It specifies which type of atom sits on each lattice site. 

• v denotes the displacement of each atom away from its ideal lattice site. 

• e is a particular electronic state (both the spatial wavefunction and spin state) when the nuclei are constrained 
to be in a state described by v. 

• E(L, cr, v, e) is the energy of the alloy in a state characterized by L, cr, v and e. 

Each summation is taken over the states that are included in the set of states defined by the "coarser" levels in the 
hierarchy of states. For instance, the sum over displacements v includes all displacements such that the atoms remain 
close to the undistorted configuration a on lattice L. 

While Equation (0) is in principle exact, practical first-principles calculations of phase diagrams typically rely 
on various simplifying assumptions. The sum over electronic states is often reduced to a single term, namely, the 
electronic ground state. The validity of this approximation can be assessed by ensuring that different structures have 
a similar electronic density of states in the vicinity of the Fermi level.0 If needed, the contribution of electronic 
entropy is, at least in its one-electron approximation, relatively simple to include without prohibitive computational 



requirements (Wolvcrton and Zungci, 1995) 



A simplifying assumption that is much more difficult to relax is the reduction of the sum over displacements v to 
a single term. This simplification has been extensively used in alloy theory, because calculating the summation over 
v involves intensive calculations. The particular displacement representing a given configuration a is typically chosen 
to be a local minimum in energy that is close to the undistorted ideal structure where atoms lie exactly at their ideal 
lattice sites. Usually, this state is found by placing the atoms at their ideal lattice positions and relaxing the system 
until a local minimum of the energy is obtained. In this fashion, the state chosen is the most probable one in the 
neighborhood of phase space associated with configuration a. In this approximation, the partition function takes the 
form of an ising model partition function: 

Z = ^^cxpHF(L,,)) (8) 

L aeL 

with E*(L, a) — min„ ie {E(L, cr, v, e)}. 

It turns out that the same statistical mechanics techniques developed in the context of the Ising model can also be 
used in the more general setting where atoms are allowed to vibrate (and where electrons are allowed to be excited 
above their ground state). All is needed is to replace the energy E*(L, a) by the constrained free energy F(L, cr, T), 
defined as: 

F(L,cr,r) = -fc s Tln(^^cxp[-^(L,cr, W , e )]) . (9) 

In other words, it is the free energy of the alloy, when its state in phase space is constrained to remain in the 
neighborhood of the ideal configuration a. This process, called the "coarse graining" of the partition function, is 
naturally interpreted as integrating out the "fa st" degrees of freedom (e.g. vibrations) before considering "slower" 
ones (e.g. configurational changes) (Cedei, 1993). This process is illustrated in Fig. |l|. The quantity to be represented 



by a cluster expansion is now the constrained free energy F(L, a, T). The only minor complication is that the effective 
cluster interactions become temperature dependent. 

There is some level of arbitrariness in the precise definition of the set of displacement v over which the summation is 
taken in Equation^. However, in the common case where there is a local energy minimum in the neighborhood of cr and 
where the system spends most of its time visiting a neighborhood that can be approximated by a harmonic potential 
well, the set of displacements over which the summation is taken has little effect on the calculated thermodynamic 
properties. Under the above assumptions, calculating the partition function of a constrained harmonic system and a 
harmonic system that allows infinite displacements gives essentially the same result: 



Un less strong electron correlation effects are present, such as charge ordering or metal-insulator transitions ([mada et al 



1998 ) , the electronic free energy can be calculated from the single electron DOS in the neighborhood of the Fermi level. 
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2 exp [-0E(L, a,v,T)]ttJ2 exp [-0E H {L,a,v, T)] 

«^ exp (L,a,u,T)] 

all n 

where E(L,a,v,T) = — ^BTln^ eglJ exp [— {3E(L, a, v, e)] and Ejj(L,a,v,T) denotes a harmonic approximation to 
E(L,a,v). In this framework, all that is needed to account for lattice vibrations, is the determination of the free 
energy of a harmonic solid in the neighborhood of any configurations a. (Appendix [f] discusses the case where the 
above assumptions are violated, that is, when no local minimum exists in the phase space neighborhood of a.) 

The problem of calculating F(L,a,T) for a set of configurations a is much more demanding than calculating the 
energy E*(L, a) for a set of a. Devising an efficient way to calculate F(L, er, T) is the fundamental problem that needs 
to be resolved in order to include vibrational effects in phase diagram calculations. 



D. Conclusion 



After presenting the basic framework enabling the calculation of the configurational free energy, this section has 
presented two important aspects of the thermodynamics of lattice vibrations. 

1. Vibrational entropy differences between phases introduce corrections to transition temperatures calculated with 
only configurational entropy. 

2. The basic alloy theory framework can be adapted to account for lattice vibrations by replacing the energy 
E*(L,a) associated with each configuration a by the free energy F(L,a,T) of a system constrained to remain 
in the phase space neighborhood of the ideal configuration a. 



III. EVIDENCE OF VIBRATIONAL EFFECTS 



In light of the large computational requirements associated with the inclusion of lattice vibrations, is it important 
to ensure that such an endeavor is worth the effort. This section reviews the experimental and theoretical evidence 
that supports the view that vibrational effects are important in the context of phase diagram calculations. There 
exists a large l iterature aimed at determining the vibra t ional properties of solids (see, for instance, (Ashcroft and 
Mcrmin , 1976 ; Maradudin et al. , 1971 ; Born and Huang , 1956| )). Here, the focus is on investigations directly related 
to the determination of vibrational entropy (or free energy) differences between phases which differ solely by the 
ordering of the chemical species on an otherwise identical parent lattice. 

This relatively narrow choice is driven by two observations. First, while there have been numerous investigations 
of the absolute vibrational properties of solids, the more difficult issue that needs to be addressed in the context 
of phase stability is the determination of accurate differences in vibrational properties. Second, it has alread y been 
estab l ished that many structural phase transformations (e.g., from fee to bec) are driven by lattice vibrations (Petry 



ct al. , 1991 ; Grimvall and Ebbsjo , 1975 ). This fact does not pose major difficulties for the purpose of phase diagram 



calculations: one can easily compute the vibrational properties of a few lattice types. The real difficulty is to calculate 
the vibrational entropy of many configurations on each of these lattices, a task which is only needed if vibrational 
properties differ substantially across distinct configurations on an identical parent lattice. 

The presentation will be mainly chronological, although deviations from that intention will be made for the sake of 
clarity. We leave a more precise description of the methods used for subsequent sections,focusing here on the results 
obtained. The key theoretical and experimental results are summarized in Tables | and III, respectively. 



A. Calculations and predicted effects on phase stability 



The idea that the state of order of an alloy could be coupled with its lattice dynamics is not new. During the 1960's, 
as the foundations of alloy thermodynamics were being established, the q uestion of the effect of lattic e vibrations was 
already being raised. Studies on the order-disorder transition of /3-brass (Booth and Rowlinson, 1955; Wojtowicz and 
1960| ), for instance, have indicated that lattice vibrations are crucial to accurately model the magnitude 



Kirkwood 



of the experimentally observed discontinuity in heat capacity at the phase transition, which determines the change in 
vibrational entropy upon disordering. 



G 



After these initial investigations, increasingly accurate models for the coupling between lat t ice v i bration a nd the 



state of order of an alloy, were then developed ( 


Moraitis and Gauticr, 


1977; 


Matthew et al. 
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1983; 


Bakker 


Bakker , 


19825; 


Bakker and Tuijn 


. L98I : 


Tuijn and Bakkei , 


1989 


; Garbulsky and Cedei , 


1994 


) . These models g 



involved unknown parameters that need to be estimated from available experimental thermodynamic data. A recurring 
theme among these studies is the idea that, for sensible choices of the stiffness of the springs connecting the atoms, 
the effect of lattice vibrations is likely to be i mportant. The estimated v ibrational entropies of diso rdering lie between 



0.05fc B, for the most conservative estimates (Moraitis and Gautier, 1977), up to the order of 0.5fcs. (Tuijn and Bakker 



1989). 



With the increased availability of computing power, the application of first-principles methods became a practical 
possibility and the unknown parameters of the theoretical models of lattice vibration have become directly com- 
putable, without relying on experimental input. Initially, only simple bulk properties, such as the bulk modulus 
were computable at a reasonable cost. This prompted the development of methods to infer vibrational properties 
f rom the knowl e dge o f elastic constants. A particularly popular sche me, the Moruzzi-Janak Schwarz (MJS) method 

Tseng 
Mohri 



( Moruzzi ct al. , 1988), was used in ma n y phase diagram calculations flAsta ct al. , 1993); Sanchez ct al. , 1991 



and Stark| |1994|; |Sanchcz and Beckcr|, |1994|; |Oh ct al.|, [1994|; |Colinct ct al.|, |1994| ; |Abrikosov et al.|, |199E 



et a|.|, |1993|;lMohri|, |1994|; |Nakamura and Mohrij |1993|; [Becker and Sanchez) , |1993| ). In the Cd-Mg flAsta et al] , [1993D , 
Ag-Cu ( Sanchez et al.| , |1991|) and Au-Ni QColinct ct al.| , |1994| j" systems, agreement with the experimentally measured 
phase diagrams was substantially improved by including vibrations in this way. In retrospect, such an improvement 
was to be expected because first-principles phase diagram calculations often greatly overestimate transitions tem- 
perature and any downward correction to the calculated transition temperatures yields improved agreement. The 
MJS approximation nearly always yields a downward correction when ordered compounds are stiffer (softer) than the 
elements in ordering (phase separating) systems, a very likely situation. 

More recently, other techniques were used to obtain vibrational properties from elastic constants. The so-called 
virtual crystal approximation (desc ribed in Appendix [C]) was use d to calculate the vibrational free energy of a dis- 
ordered alloy in the Ni-Cr system (Craievich and Sanchez. 1997). The calculated vibrational free energy exhibited 
qualitatively the same concentration-dependence as the vibrational free energy obtained by subtracting the experi- 
mentally determined free energy from the calculated configurational free energy. The virtu al crystal approximation s 
was also used in calculations of the lattice dynamics of ordered and disordered C U3AU (Cleri and Rosatc , 1993 ). 
These calculations, which relied on a tight-binding Hamiltonian (see, for instance, ( Harrisorj ~1989; Pettifoi, 1992)), 
predicted a 0.12kg increase in the vibrational entropy upon disordering. 

In view of the large computational requirements of accurate ab initio methods, many researchers have sought to 
calculate vibrational properties with simpler energy models, whose lower computational requirements enable a more 
accurate handling of issues such as an harmonicity or the repr esentation of the disordered state. The development of 
the embedded atom method (EAM) ( Daw and Baskes , 1984 ) offered the opportunity to accurately model metallic 
alloys at a reasonable computational cost. Investigations based on the EAM have typically found a large vibrational 
entropy change upon d isordering in metallic alloys. The vibrational entropy change for Cu.^ Au was p r edicted to be 
O.lOfc a (|Ackland , 1994), w hile for M3AI, values ranging from 0.22fcs to 0.29fcs were obtained ( lAcklancj |l994| ; |Althohj 
et a[j, |l997j ; [Althoff et al.| , |l998j ; [Ravelo et al] , |1998j )p|. 

Other researchers c onstructed pair potentials from an equation of state determined from first-principles calculations 
( |5haojun et~aT[ , |1998[ ). The resulting pair potentials were then used to calculate disordering vibrational entropies with 



no further approximations. This method attributed a vibrational entropy change of 0.11 kg to the disordering reaction 
of the FeaAl compound. 

Rather unexpected results were uncovered as it became possible to compute vibrational entropy differences from a 
complete la t tice d ynamics analysis as well as state-of-the-art ab initio techniques. Calculations on the Si-Ge system 
(Garbulsky, 1996| ) found almost no effect of lattice vibrations: the vibrational entropy of formation of the metastable 
zincblende structure was a mere — 0.02fcs- The first ab initio calculation of a vibrational entropy of disordering 



(van de Walle et al. , 1998) placed an upper bound of 0.05fcs i n the case of the order-disorder transition of the N13AI 



compound, in sharp contrast with previous EAM calculations ( Ackland , 1994 ; Althoff ct al. , 1997 ; Althoff ct al. , 1998) 
which found a much larger value. Although the disagreement simply originated from a difference in the predicted 
volume expansion of the alloy upon disordering, the difference between EAM and ab initio results docs indicate 
that vibrational entropy differences are very sensitive to the e nergy model used. In another int crmetallic compound, 
Pd 3 V, the vibrational entropy change upon disordering of the (van de Walle and G.Ceder, 2000) was calculated to be 



5 Disordered M3AI is actually a metastable phase. The values quoted here are at the highest temperatures reported by the 
investigators, as close as possible to the true disordering temperature that would be observed if the alloy did not melt before. 
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— 0.07/cb, although the simplest theories put forward in the earliest investigations of lattice vibrations in alloys would 

predict this value to be large and positive. 

In the first phase diagram cal culation based on a full lattice dynamics analysis (Tepesch et al., 1996| ), the SCPIB 



method ( Boyer and Mchl , 1993 ) was used to calculate that vibrational effects lower the top of the miscibi lity gap 
by abo ut 10% in the CaO-MgO system. Subsequently, first-principles calculations on the Cu-Au system (Ozolini 



et aQ |!998cD reported a reduction of a about 20% in the transition temperatures. However, in both cases, the 
resulting correction on the phase boundaries decreased the agreement with the experimentally determined phase 
diagram, suggesting that other potential sources of error, such as the precision of the energy method used or the 
convergence of the cluster expansion, have to be investigated. 

Very recently, two rather striking examples of systems where lattice vibrations do have an important effect and 
where their inclusion results in a a dramatically improved agreement with experimental measurements have been 



uncovered. Lattice vibrations are responsible for a 27- fold increase in the solubility of scandium in aluminum ( Ozolini 
and Asta, [2001 ), resulting in a nearly perfect agreement with e xperimentally determined solubility limits. Vibrational 
contributions were also shown (Wolverton and Ozolins, 2001) to be essential to correctly predict the relative stabil- 
ity of the different precipitates which constitute the well-known Guinicr-Prcston zones responsible for precipitation 
hardening in the Al-Cu system. It is interesting to note that these two successful examples did not require the use of 
a cluster expansion and did not necessitate a large number of separate first-principles calculations. Each calculation 
could therefore be carried out with an extremely high precision. When all other sources of errors are well controlled, 
the effect of lattice vibrations can be more accurately quantified. These examples offer good hope that ab-initio 
calculations of a complete phase diagram that include lattice vibrations can provide an accuracy comparable to ex- 
periments, provided that a well converged cluster expansion can be obtained and that highly accurate first-principles 
calculations are used to construct it. 



B. Comparison with Experiments 



Over the last 10 years, advances in experimental techniques have made it possible to directly measure vibrational 
entropy differences, instead of inferring them from discrepancies between measured thermodynamical data and cal- 
culated estimates of configurational contributions to the free energy. Experimental investigations have thus provided 
independent assessments of the role of lattice vibration on phase stability. 

The first dire c t mea surement was obtained from differential calorimetry measurements on the N13AI compound 
( Anthony et al. , 1993 ) and found the vibrational entropy of disordering to be at least 0.19/ cb- This finding was 



corroborated with subsequent incoherent neutron scattering measurements (Fultz et al., 1995), which bracketed its 



value between O.l fcg and . 3fcg. These findings fuel e d much of the inte r est of the recent theoretical literature on th e 
Ni 3 Al compound ( |Ackland| , |1994| ; [Althoff et al] , [1997| ; [Althoff et al] , |1998j ; |Ravelo et al] , |1998| ; [van dc Wallc et al] , [1998p . 
Unfortunately, the agreement between experimental and theoretical determinations is relatively poor. Even among 
studies in which the magnitude of the vibrational entropy is similar, its proposed physical origin differs substantially: 
according to experiments, t he vibrational entropy change occured with es sentially no change in volume, while most 



calculations ( Althoff et al. , 1997 ; Althoff et al , 1998 ; Ravelo et al. , 1998j ) attribute the vibrational entropy change 



almost entirely to a volum e change. Many of these conflicting findings were clarified by first-principles calculations 



( van de Walle et al. , 1998 ), which predicted both a very small volume change and a very small vibrational entropy 
change upon disordering. These results indicated that (i) the apparent large vibrational entropy change observed 
experimentally could very well be entirely explained by the nanocrystalline nature of the samples and the use of the 
virtual crystal approximation in interpreting the data and (ii) the large volume expansion upon disordering predicted 
by the EAM is probably an artifact of the method, as it would have otherwise been clearly visible experimentally. 
To be fair, Fultz et al. and Althoff et al. were fully aware of these potential problems; Ab initio calculations merely 
made it possible to better quantify the errors introduced. The general consensus among researchers is now that given 
the numerous difficulties faced when studying N13AI, unambiguous evidence of the importance of lattice vibrations 
should probably be sought in other systems. 



Sim ilar c alorimetry measurements were then performed on the FesAl ( Anthony et al. , 1994| ) and CU3AU (Nagcl 



ct al., 1995) compounds, where more conclusive results could be obtained. The vibrational entropy of disorderin g of 
FeaAl was determined to be O.lOfcs, a result which was later corroborated by calculations (Shaojun ct al. , 1998). In 



the case of CU3AI the experimental resu l t, (0.1 4 ± 0.05) ks, showed very good agreement with the earli e r theo retical 



predictions of 0.12/cb ( Cleri and Rosatc] , 1993 ). Subsequent linear-response calculations ( Ozolins et al. , 1998c ) yield 



values ranging from 0.06 to 0.08fce (depending on temperature), which also compares favorably with the experimental 
results. 

The estimation of vibrational effects was also addressed by directly probing the lattice dynamics through neutron 
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scattering measurements. The vibrational entropy of formatio n of a disordered a lloy in the Fe-Cr system was obtained 
from single crystal measurements of phonon dispersion curves ( Fultz ct al. , 1995| ) in the virtual crystal approximation. 
Values ranging from 0.14fcs to 0.21fce were obtained, depending on concentration. 

In order to determine the lattice dynamics of disordered alloys beyond the virtual crystal approximation, the 



inco heren t neutron scattering technique was extens i vely used and refi n ed ( Nagel et al 
nfl|997 



so 



Robertson et al 



1996 



Nagel, Fultz and Robert 



1999J; [Bogdanoff et al] , |l999| |Fultz et al] , |1995| ). With this technique, the analysis of the 



experimental data is considerably simplified when the species present have comparable incoherent neutron scattering 
intensities, which lead to the study of two compounds satisfying this requirement: N13V and C03V. Measurements 



on the M3V compound (Nagel et al., 1996) found a surprisingly small vibrati onal entropy change upon disorder ing 



0.04fcs, while the C03V compounds exhibited a relatively large value 0.15fcs ( Nagel, Fultz and Robertson , 1997] ). A 
related stu dy found the vibration al entropy change associated with the fee- hep transition of the C03V compound to 
be 0.07/cb (Robertson et al., 1999). It is interesting to note that the disordering reaction exhibits a larger vibrational 



entropy change than the allotropic transformation in C03V. Perhaps more importantly, the investigations of the M3V 
and C03V compounds presented the first experimental evidence of important anharmonic effects. 



The same technique of incoherent neutron scattering was employed to revisit the Cu-Au system ( Bogdanoff ct al 



1999). The vibrational entropy of formation of the CU3AU compound w as found to be 0.06^5 at 300 K, corroborating 

At 800^, the 



1999|) . 



earlier estimations based on phonon dispersion curve measurements ( Bogdanoff and Fultz . 
measu red va lue of (0.12 ± 0.04)/cb is consistent with the value of 0.20/cs obtained with ab-initio calculations (Ozolini 
ct aQ, |l998cj) . 

The vibrational entropy of formation of v arious ordered compounds, obtained from single crystal phonon disper- 
sion measurements, were recently compiled (Bogdanoff and Fultz, 1999) and show formation values of up to 0.5/c£. 
However, this compilation contains many systems where the alloy has a crystal structure that differs from the one 
of the pure elements and the formation values thus also include the vibrational entropy change associated with a 
structural transition. When all these cases are excluded, the maximum vibrational entropy change decreases to a 
more conservative upper bound of 0.20fcs, which is reached in the ordered phase of M3AI. 



C. Conclusion 



Although early investigations of the impact of vibrational effects on phase stability consistently found large effects, it 
is now becoming apparent, as more precise theoretical and experimental techniques became available, that vibrational 
effects are often, but not always, large. It is therefore important to identify the factors which determine when they 
are, so that the effort devoted to calculating them is proportional to their expected magnitude. 



IV. THE ORIGIN OF VIBRATIONAL ENTROPY DIFFERENCES BETWEEN PHASES 



We have presented the framework that allows for the inclusion of vibrational effects in phase diagram calculations. 
However, the formalism presented so far does not directly provide any intuition regarding the origin of vibrational 
entropy differences. This intuition is important to be able the predict when vibrational effects should be important 
and, when they are, which approximation should be used to calculate them. 

Three mechanisms have been suggested to explain the origin of vibrational entropy differences in alloys. We will 
discuss them in turn. 



A. The "bond proportion" effect 



In most theoretical studies based on simple models systems ( 


Vlatthew et al. 




1983; Bakker 




Tuijn 




1986; Tuijn and Bakker, 


1989; 


Garbulsky and Cede: 


, 1994 


), the effect of the state of orde 



1982q; Bakker and 



vibrational entropy has been attributed to the fact that bonds between different chemical species have a different 
stiffness than the bonds between identical species. When the proportion of each type of bond in the alloy changes 
during, for instance, an order-disorder transition, the average stiffness of the alloy changes as well, resulting in a 
change of its vibrational entropy. This so-called "bond-proportion" mechanism is illustrated in Fig. ^ in the case 
of an order-disorder transition. In a system with ordering tendencies, the bond between unlike atoms are associated 
with an increased stability and are thus expected to be stiffer than bonds between alike atoms. Since disordering 
reduces the number of bonds between unlike atoms in favor of bonds between similar atoms, the disordered state 
is expected to be softer, and thus, have a large vibrational entropy. Vibrations would then tend to stabilize the 
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disordered state relative to the ordered state, reducing the transition temperature. A similar reasoning in the case of 
a phase separating system shows that the disordered state should be softer than a phase separated mixture, indicating 
that the miscibility gap should be lowered as a result of vibrational effects. 

The presence of a "bond proportion" effect can be readily identified from the nature of the changes taking place in 
the phonon densities of states during an order-disorder transition. In the ordered alloy, the very stiff nearest-neighbor 
bonds should be associated with high frequency optical modes peaks. As the alloy disorders, the height of these peaks 
should decrease since the number of stiff bonds decreases. This characteris tic signature of the "bond proportion" 



mechanism in the phonon DOS has bee n repeatedly o b serve d in experiments ( Fultz et al. , 1995 ; Anthony et al. , 1994 ) 
as well as in theoretical calculations ( Althoff et alj , 1997) (although, in the latter study the "bond proportion" 
mechanism was dominated by other effects discussed below). 



B. The volume effect 



It is well known that vibrational entropy of a given compound varies with volume: This dependence is responsible 
for thermal expansion. It is thus expected that the volume change that typically occurs during solid-state transitions 
should also result in a change in vibrational entropy. As the alloy expands (or contracts), as a result of a change in its 
state of order, the stiffness of all chemical bonds decreases (or increases) . The resulting change in vibrational entropy 
is entirely due to anharmonicity, in contrast to the "bond proportion" effect. When the volume mechanism operates 
alone, the phonon DOS should exhibit an overall shift when the volume changes (See Fig. ||). 

This shift is usually accompanied by a change in the shape of the phonon DOS, so that visual inspection of the 
phonon DOS is usually not sufficient to identify the shift due to the effect of volume. Theoretical investigations of 
this effect have thus relied on a simple thought experiment consisting in separating the vibrational entropy change 
upon disordering at constant volume from the vibrational entropy change resulting solely from the volume expansion 
of the disordered state. In this fashion, the importance of volume changes was first obs erved in embed ded atom 
method (EAM) calculations of disordering reaction of the M3AI and CU3 AU compound s ( Ackland , 1994 ) an d late r 



corroborated by subsequent EAM calculations on the Ni 3 Al compound (Althoff et al., 1997; Ravelo et al. , 1998) 



These more recent calculations also found that the volume effect is magnified by the fact that the linear thermal 
expansion coefficient of different phases can differ substantially (by about 5 x 10~ 6 K^ 1 relative to an absolute value 
of about 15 x 10 -6 K^ 1 ). First-principles calculations on the Cu-Au system (Ozolins et al. , 1998c ) revealed a similar 
finding. I n con trast, first-principles calculations on the M3AI (van de Walle et al., 1998) and PdsV (van dc Walle and 
G.Ceder , 2000 ) compounds found only a small difference between the thermal expansion coefficients of the ordered 
and the disordered phases (less than 1 x 10~ 6 K^ 1 ). 

Even though the temperature-dependence of vibrational entropy can be large, it is important to keep in mind that 
a given vibrational entropy difference arising from anharmonic contributions has a smaller impact on the vibrational 
free energy than a harmonic contribution of the same magnitude. The reason is that anharmonic contributions 
to the vibrational entropy are always partly canceled by anharmonic contributions to the vibrational enthalpy. In 
contrast, such cancellation does not occur for harmonic contributions.^ The quasiharmonic approximation pred icts 



19984) 



that anharmonic contributions have exactly half the effect of harmonic contributions (see ( Ozolins et al. 
Appendix ^) when the volume dependence of the energy is quadratic while the volume dependence of the vibrational 
entropy is linear. 

Experimental measurements have not identified the volume mechanism as a major source of entropy differences 
since the simple thought experiment that allows its identification in calculations cannot be performed experimentally. 



The effect of thermal e xpan sion on the phonon DOS, however, is clearly seen experimentally (Nagel et al., 1996; Nagel 



Fultz and Robertson, 1997), due to the fact that thermal expansion causes shifts in the phonon DOS that arc not 



accompanied by subs tantial chang es in its shape. Measu rements on M3V (Nagel et al., 1996), C03V (Nagel, Fultz 



and Robertson, 1997) and CU3AU (Bogdanoff et al., 1999| ) all show that anharmonic contributions are not negligible 



6 A temperature-dependence of vibrational entropy necessarily introduces a temperature-dependence of the vibrational en- 
thalpy, as a consequence of the following thermodynamic relation: dI ^ ib — T d g^i b . The anharmonic vibrational free energy is 
then a sum of two competing contributions: F v y, = H V ib — TS V ib- In contrast, harmonic contributions to the vibrational en- 
thalpy are configuration-independent in the high-temperature limit, by the equipartition theorem, and give no net contribution 
to vibrational free energy differences. Vibrational entropy differences originating from harmonic contributions thus enter the 
expression for vibrational free energy differences directly, without any partial cancellation from the enthalpic term. 
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C. The size mismatch effect 



The third advocated source of vibrational entropy changes is the effect of atomic size mismatch. When atoms of 
different sizes are constrained to coexist on a lattice, the atoms can experience compressive (or tensile) stress that 
results in locally stiffer (or softer) regions. When large atoms sit on neighboring lattice sites, the amplitude of their 
vibrations is reduced, i.e., the alloy tends to be locally stiffer. Conversely, when small atoms sit on neighboring lattice 
sites, the extra room available results in a locally softer region. 



The phenomenon was first noted in EAM calculations on the CU3AU compound (Ackland, 1994), where, in the 



disordered state, the presence of highly compressed pairs of Au atoms lowers the vibrational entropy of the disordered 



state. A similar effect was found in first-principles calculations on M3AI (van de Walle et al. , 1998), whe re very 
compressed pairs of Al ato ms where found in the disordered state. In first-principles calculations on Pd3V (van dc 
Walle and G.Cedei , 2000 ), an even more intriguing size-related effect was observed: all three types of chemical 
bonds have incompatible equilibrium lengths and the vibrational entropy changes can be explained solely by the large 
relaxations of th e atoms away from t heir ideal lattice sites in the disordered state. A summary of these observations 
can be found in (Morgan et al., 2000] ). 



In the experimental literature, the size mismatch effect has been described with the help of a stiff sphere picture first 
introduced in an investigation of the CU3AU compound (Nagel et al., 1995| ). The fundamental intuition behind this 
picture is that a chemical bond becomes stiffer when the two bonded atoms have touching atomic "spheres" . Further 
evidence for this stiff sph ere picture was provided by a systematic analysis of the vibrational entropy of formation of 
various LI2 compounds ( Bogdanoff and Fultz , 1999 ), which shows a correlation with the difference in radii between 
the two alloyed species. 



V. COMPUTATIONAL TECHNIQUES 



To understand the nature of the difficulties encountered, it is instructive to first consider how, in principle, the 
vibrational properties of a single configuration a can be calculated with an arbitrary accuracy. The techniques 
present ed in this section are the tools that were used to investigate the importance of lattice vibrations presented in 
Section [ill A . 

Phase diagram calculation involves computing vibrational properties for a set of configurations a. Carrying out 
the full phonon problem for each configuration results in undue computational requirements. Nevertheless the formal 
solutions presented here play an important role in devising practical ways to include vibrational effects in phase 
diagram calculations. This section first focuses on the treatment lattice vibration within the harmonic approximation, 
before addressing the issue of anharmonicity. Finally, important consideration concerning the energy models used as 
an input for these procedures are discussed. 



A. Lattice vibrations in the harmonic approximation 



In this section, we review the problem of determining the constrained free energy of a system in the neighborhood 
of a configuration er, under the assumption that the system spends most of its time in a region near a local energy 
minimum, where a harmonic approximation to the energy surface is accurate . In th is approximation, the fr e e ene rgy 
determination reduces to the well-known phonon problem (Maradudin et al., 1971), (Ashcroft and Mcrmin, 1976). 



1. Theory 



Consider a system consisting of N atoms. Let Mi be the mass of atom i and u (i) be its displacement away from its 
equilibrium positions. Time derivatives are denoted by dots while Greek letter subscripts denote one of the cartesian 
components of a vector. In the harmonic approximation, the energy of the system can be written as: 



(10) 



where 



d 2 E 



du a (i) dufj (j) 



(11) 



u(l)=0 VI 
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The 3x3 matrices <E> are called the force constants tensors, as they relate the displacement of atom j to the 
force / exerted on atom i: 



f(i) = $(i,j)u{j) 



(12) 



Such a harmonic approximation of a solid is often referred to as a Born-von Karman model. 

Note that contrary to usual treatment, we do not immediately impose translational symmetry, in order to derive a 
few general results that also apply to systems such as disordered alloys. 

The substitution e (i) = \fWlu (i) yields: 



1 



1,3 



■U) 



(13) 



The 3iV eigenvalues A m of the matrix^ 



/ JHi 



D 



\ 



<S>(N,1) 
VA/ivA/i 



•S>(1.N) 

VMiMjv 



<t>(N,N) 
\/M N M N 



(14) 



then give the frequencies v m = \/ X m of the normal modes of oscillation. In the harmonic approximation, the 
knowledge of these frequencies is sufficient to determine the thermodynamic quantities we are interested in. This 
information is conveniently summarized by the phonon density of states (DOS), which gives the number of modes of 
oscillation having a frequency lying in the interval [v, v + dv\ : 



1 3N 



v m ) ■ 



(15) 



I t can be shown that th e free energy of the system (restricted to remain close to a given configuration a) is given by 
( [Maradudin ct aD , |l97l| ): 



F E* f°°, ( , / hv 



N N 
E* 
~N 



N 



^2 m ( 2 smn 



2k B T 



g (v) dv 



where E* is the potential energy of the system at its equilibrium position and h is Planck's constant. Phase tran- 
sitions in alloys typically occur at a temperature where the high temperature limit of this expression is an accurate 
approximation: 



F E* , 

— = h k B T 

N N 



In 



h v 

( hV; 



N N n yk B T 



g (u) dv 



The usual criterion used to determine the temperature range where high temperature limit is reached is the Debye 
temperature. Note that the factor -^-^ is often omitted because it cancels out when calculating vibrational free energy 
differences. In the high temperature limit, another important f o rm of cancellation occurs: The at omic masses have 
no effect on the free energies of formation (Grimvall and Rosen, 1983; Garbulsky and Ceder, 1994). This important 



7 Among the 3iV eigenvalues, the 6 eigenvalues associated with rigid body translations and rotations are zero. In the thermo- 
dynamic limit, these few degrees of freedom are inconsequential. To avoid notational complications, we simply assume that the 
solid is fixed to a reference frame by springs so that the resulting dynamical matrix has 37V nonzero eigenvalues. 
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result, shown in Appendix [A|, rules out that masses play any significant role in determining phase stability at high 
temperatures. 

As mentioned before, a convenient measure of the magnitude of the effect of lattice vibrations on phase stability is 
the vibrational entropy, which can be obtained from the vibrational free energy by the well known thermodynamical 
relationship S v q, = — ^dT* . Contrary to the vibrational free energy of formation, the vibrational entropy of formation^ 
is temperature-independent in the high-temperature limit of the harmonic approximation, allowing a unique number 
to be reported as a measure of the importance of vibrational effects. 

In a crystal, the determination of the normal modes is somewhat simplified by the translational symmetry of the 
system. Let n denote the number of atoms per unit cell. Let u (') denote the displacements away from its equilibrium 

position of atom i in cell I. Let ^) be the force constant relative to atom i in cell I and atom j in cell I' and let 

Bloch's theorem indicates that the eigenvectors of the dynamical matrix are of the form 



e z.27r(fc-Z) e 



(16) 



where / denotes the cartesian coordinates of one corner of cell I and k is a point in the first Brillouin zone. This 
fact reduces the problem of diagonalizing the 3N x 3N matrix D to the problem of diagonalizing a 3n x 3n matrix 
D (k) for various values of k. This can be shown by a simple substitution of Equation ( |l6| ) into Equation (13)|] The 
dynamical matrix D (k) to be diagonalized is given byp*| 



/ *(?0 



D(k) = J2^ Mk ' l) 
I 



v/A/iMi 



, *(nl) 

\ \/M„ jl/i 



HZ) 

s/M n M n 



(17) 



/ 



As before, the resulting eigenvalues (fc) for % — 1 . . . n, give the frequencies of the normal modes {vi (k) =-^\J\ (&))• 
The function (k) for a given i is called a phonon branch, while the plot of the fc-dependence of all branches along 
a given direction in k space is called the phonon dispersion curve. In periodic systems, the phonon DOS is defined as 



3n „ 

(v) = J2 / 8(v- Vi (k)) dk 

J BZ 



(18) 



where the integral is taken over the first Brillouin zone. 



2. Force Constant Determination 



The above theory relies, of course, on the availability of the force constant tensors. The determination of these 
force constant tensors is the focus of this section. Before describing the methods used for their determination, we will 
first review important properties of the force constant tensors. 

While the number of unknown force constants to be determined is in principle infinite, it can, in practice, be 
reduced to a manageable finite number with the help of the following two observations. First, the force constant 
$(i,j) between two atoms i and j beyond a given distance can be neglected. Second, the symmetry of the crystal 
imposes linear constraints between the elements of the force constant tensors. 

The accuracy of the approximation made by truncating the range of force constant can be tested by gradually 
increasing the range of interactions, until the quantities to be determined no longer vary substantially. It is important 
to note that most thermodynamic quantities can be written as a weighted integral of the phonon DOS and their 



The absolute value of the vibrational entropy is not constant at high temperature, but its temperature-dependence does not 
vary across distinct phases and thus formation values are temperature-independent. 
9 And changing the summation over atoms by summations over atoms and cells. 

10 The reader should be aware that there are many possible conventions regarding the phase factor: for instance, e* 2 "''' 1 ', 
e -i2n(k-i)^ e »(fc-O j e '2n(k-(i+x(j))) ^ere x(j) is the coordinate of atom j within the cell. While all convention yield different 
dynamical matrices, they all have the same eigenvalues. 
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convergence rates are thus muc h faster than the pointwise convergence rate of the phonon DOS itself (Garbulsky 



1996 ; van de Walle et al. , 1998 ). That is, the errors on the DOS at each frequency tend to be quickly averaged out 



when the contributions of each frequency are added. 

The restrictions on the force constants imposed by the symmetry of the lattice can be expressed as follows. Consider 
the force constant of atoms i and j located at x(i) and x(j) and consider a symmetry transformation that 

maps a point of coordinate x to Sx + 1, where S is a 3 x 3 matrix and t and 3x1 translation vector. In general, if 
the crystal is left unchanged by such a symmetry operation, the force constant tensors should be left unchanged as 
well. This fact imposes the following constraints on the spring tensors: 



Sx(i)+t = x(i') 
Sx(j)+t = x(j>) 



<S>(i,j) = S T <S>(i',j')S. 



(19) 



Additional constraints on the force constants can be derived from simple invariance arguments. The most important 
constraints, obtained by noting that rigid translations and rotations must leave the forces exerted on the atoms 
unchanged, are 

$(»,«) = - 



Additional constrains can be found in (Maradudin ct al., 1971; Born and Huang, 1956| ) 



There are essentially three approaches to determining the force constants: analytic calculations, supercell calcula- 
tions and linear response calculations. Analytic calculations are only possible when the energy model is sufficiently 
simple to allow a direct calculation of the second derivatives of the energy with respect to atomic displacements, as in 
the case of empirical pair potential models. For first-principles calculations, either one of the two following methods 

have to be used. 

a. The supercell method The supercell method ( Wei and Chou| , 1992), (Garbulsky, 1996| ) consists of slightly per- 
turbing the positions of the atoms away from their equilibrium position and calculating the reaction forces. Equating 
the calculated forces to the forces predicted from the harmonic model yields a set of linear constraints that allows the 
unknown force constants to be determined .p] 



F(i) = $(i,j)u(j) 



(20) 



When the force constants considered have a range that exceeds the extent of the primitive cell, a supercell of the 
primitive cell has to be used. (The simultaneous movement of the image atoms introduces linear constrains among 
the forces that prevent the determination of some of the force constants.) 

While any choice of the perturbations that allows the force constants to be determined is in principle equally valid, 
a few simple principles drastically narrow down the number of perturbations that need to be considered. For a given 
supercell, there is a only of finite number of non-redundant perturbations to consider. 

A minimal set of non-redundant perturbations can be obtained as follows. 

• Consider in turn each atom in the asymmetric unit of the primitive cell. 

• Mark the chosen atom (and its periodic images in the other supercells) and consider it as distinct from other 
atoms of the same type. (This operation effectively removes some of the symmetry operation of the space group 
of the crystal.) 

• Construct the point group {Si} of the site where this atom is located. (Si is a 3 x 3 matrix.) 

• Move the chosen atom along a direction u\ such that the space spanned by the vectors SiU\ (for all i) has the 
highest dimensionality possible. 



11 Equalities between calculated and predicted energies can be used as well. Using energies alone to determine the force con- 
stants would be a rather inefficient use of the information provided by ab-initio calculations. Once a first-principles calculations 
of the energy of a distorted structure has been completed, the calculation of the forces acting on the atoms is computationally 
inexpensive. The knowledge of the energy provides a single equation while the knowledge of the forces provide up to 3 equations 
per atom. 
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• If the resulting dimensionality is less than three, consider an additional direction U2 such that the space spanned 
by the vectors SiUj for j — 1,2 has the highest dimensionality possible. 

• If the resulting dimensionality is less than three, consider a direction 1*3 orthogonal to u\ and ui- 

The resulting displacements Uj for all atoms in the asymmetric unit gives a minimal list of perturbations that is 
sufficient to find all the force constants that can possibly be determined with the given supercell. This result follows 
from the observation that any other possible displacement can be written as a linear combination of the displacements 
considered above (or displacements that are symmetrically equivalent to them). 

When determining force constants with the supercell method, it is important to verify that the presence of small 
numerical noise in the calculated forces does not result in too much error in the fitted force constants. To minimize noise 
in the fitted force constants, it may be necessary to use more than the minimum possible number of perturbations. The 
additional perturbations should ideally be based on different supercells, to minimize the systematic errors introduced 
by the movement of the image atoms. 

When ab-initio calculations are used to calculate the forces, it is especially important to iterate the electronic self- 
consistency steps to convergence. Even though the energy may appear to be well converged, the forces may not yet 
be. Energy is the solution to a minimization procedure, while forces are not. As a result, errors on the energy are of a 
second order in the minimization parameters, while the errors on the forces are of the first order in the minimization 
parameters. For the same reason, special attention should be given to the structural relaxations. 

The true system is not exactly harmonic and the calculated forces may exhibit anharmonic components that 
introduce noise into the fitted force constants. This problem can be alleviated by considering an additional set of 
perturbations, where the displacements have the opposite sign. Subtracting the calculated forces obtained for this new 
set of displacements from the corresponding displacements of the opposite sign exactly cancels out all the odd-order 
anharmonic terms. Of course, for perturbations such that the negative displacement is equivalent by symmetry to 
the corresponding positive displacement, this duplication is unnecessary, because the terms of odd order are already 
zero by symmetry. 



1997: Wei and Chou, 1992: Gar- 



Additional guidelines for fitting force constants can be found in (Ackland et al 
bulsEyj |1996| ). 

b. Linear response Linear response calculations seek to directly evaluate the dynamical matrix for a set of fc 
points. The starting point of the linear response approach is evaluation of the second-order change in the electronic 
energy induced by a atomic displacements from pe rturbation theory Within this framework, practical schemes to 



compute vibratio nal properties in semi conductors (Baroni ct al. , 1987 ; Giaunozzi ct al., |l991 ; Gonzc ct al. , 1992 



Waghmarc, 1996) and metallic systems (dc Gironcoli 



199E; Quong and Klein, 1992; Ozolins, 1996) have been devised 



In this section we will not discuss the theory behind linear response calculations which can be found in a recent review 



(Gonze, 1997; Gonze and Lee, 1997), but rather focus on how the results of linear response calculations can be used 



in the context of alloy phase diagram calculations. 

The dynamical matrices calculated from linear response theory are exact in the sense that they account for arbitrarily 
long-range force constants. While in the supercell method inaccuracies arise from the truncation of the force constants, 
the limit in precision for linear response calculations arises from the use of a small number of fc points to sample the 
Brillouin zone. To address this issue, two methods can be used. 



A set of special k points c an be selected through the Chadi-Cohen ( Chadi and Cohen , 1973 ) or Monkhorst-Pack 
(Monkhorst and Pack, 1976) schemes. Special k points are selected so that the integral over the Brillouin zone of a 
function h (fc) that contains no Fourier components above a given frequency can be exactly evaluated by a weighted 
average of the function at each special point. Since thermodynamic quantities can be written as integrals of functions 
of the dynamical matrix / (D (fc)) over the Brillouin zone, the procedure is straightforward to apply in this context. 
The other app roach is the so-called Fourier inversion method (see, for instance, ( piannozzi ct alj , 1991; Quong 



and Klein, 1992j )). The calculated dynamical matrices from a set of k points are used to determine the value of the 



force constants up to a certain interaction range. The resulting harmonic model can then be used to calculate the 
dynamical matrix at any point in the Brillouin zone, allowing a much finer sampling of the Brillouin zone for the 
purpose of performing the numerical integration required to determine any thermodynamic quantity. 

The Fourier inversion method is preferable when the function / (D (fc)) to be integrated exhibits high-frequency 
components, while the dynamical matrix itself, D (fc), does not. Such a situation would arise when the function / is 
highly nonlinear. The smoothness of D (fc) then ensures that it can be represented with a small number of Fourier 
components. The less well-behaved function / (D (fc)) can then be accurately integrated with as many fc points as 
needed, using the dynamical matrix D (fc) calculated from the spring model. 

In the case of vibrational free energy calculations, the special fc points method has been observed to converge 
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rapidly with respect to the number of k points (Garbulsky and Ceder, 1996[ )p|, so that the Fourier inversion method 
is probably unnecessary.^ 

For a given set of special k points, there is an approximate correspondence between the number of Fourier compo- 
nents that can be integrated exactly and the range of force constants that can be determined. The correspondence is 
exact only when the lattice has one atom per cell and when the function / is linear 

While supercell and linear response calculations are in principle equivalent in terms of the information they provide, 
they have complementary advantages in terms of computational efficiency. The linear response method is the most 
efficient way to perform high-accuracy calculations that would otherwise be tedious and computer intensive with the 
supercell method. However, when a high accuracy is not needed, the supercell method has the advantage that various 
simplifying assumptions regarding the structure of the force constant tensors can transparently be used to drastically 
reduce computational requirements. It is not clear at which level of accuracy the cross-over between the efficiency 
of each approach occurs, but it is important to keep both approaches in mind. Another consideration is that in the 
continuously evolving field of computational solid state physics, new first-principles energy methods are continually 
developed, and the derivation of the appropriate linear response theory always follows the derivation of simple force 
calculations. Hence, despite the elegance of linear response theory, it is to be expected that the supercell method will 
always remain of interest. 



B. Anharmonicity 

While the harmonic approximation is remarkably accurate given its simplicity, it has one important limitation: It is 
unable to model thermal expansion and its impact on vibrational properties. Both the free energy F and the entropy 
S can be obtained from the the heat capacity C p : 



Q-dT and F = E\ 
T 



T=0 



SdT, 



(21) 



Hence, a simple way to account for thermal expansion is to use the following well known thermodynamic relationship 
between the heat capacity at constant pressure C p and at constant volume C v : 



C p = C v + BVTa 2 



(22) 



where a is the coefficient of volumetric thermal expansion while B is the bulk modulus. In a purely harmonic model, 
there is no thermal expansion and C v is equal to C p . The term BVTa 1 can thus be viewed as correction arising from 
anharmonic effects. 

Equatio n (|22]) is direc t ly us e ful in the context of exp erimental measurements where C p , V and a can be directly 



measured (Nagel ct al., 1996; Robertson et al., 1999| ). In the following section, we describe the computational 



techniques used to handle anharmonicity. 



1. The quasiharmonic model 



A simple modification to the harmonic approximation, called the quasiharmonic approximation, allows the cal- 
culation of thermal expansion at the expense of a moderate increase in computational cost. In the quasi-harmonic 
approximation, the phonon frequencies are allowed to be vol ume-depe n dent , which amounts to assuming that the force 
constant tensors are volume-dependent (see, for instance, (Grimvall, 1986| )). This approximation has recently been 
shown to be extremely reliable, enabling accurate first-principles calculations of the thermal expansion coefficients of 
many elements up to their melting points (Quong and Lui, 1997). 



12 Calculations of the authors, based on data from ( jOzolins et al. , 1998c) also support this finding. 

13 Note that the function whose integral gives the free energy exhibits a logarithmic singularity at the T point, which could lead 
to high frequency components that are difficult to integrate accurately. However, in three-dimensional systems, this logarithmic 
singularity contributes very little to the value of the free energy, so that the rate of convergence of the integral as a function 
of the number of k points is not dramatically slowed down by the presence of the singularity. As a result, the special k point 
method can safely be used in practical calculations of the vibrational free energy. 

14 As can be shown by a simple Fourier transform. 
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The best way to understand this approximation is to study a simple model system where it is essentially exact. 
Consider a linear chain (with periodic boundary conditions) of identical atoms interacting solely with their nearest 
neighbors through a pair potential of the form: 

U (r) = air + a 2 r 2 + a 3 r 3 . (23) 

Let L be the average distance between two nearest neighbors and let u (i) denote the displacement of atom i away 
from its equilibrium position. The total potential energy (per atom) of this system is then given by 

— = — «i ( L + u (*) " u (* + 1)) + 02 {L + u (i) - u (t + l)) 2 

i 

+a 3 (L + u(i) -u(i + l)f 

This expression can be simplified by noting that all the terms that are linear in (u (i) — u(i + 1)) cancel out when 
summed over i. 



a x L + a 2 L 2 + a 3 L 3 + (a 2 + 3a 3 L) (it («)-«(« + ^ 

i 

+0((u (i)-u{i + l)f^j 



The first three terms, a\L + a 2 L 2 + a 3 L 3 , give the elastic energy of a motionless lattice while the remaining terms 
account for lattice vibrations. The important feature of this equation is that, even within the harmonic approximation, 
the prefactor of the harmonic term, (a 2 + 3a 3 L), depends on the anharmonicity of the potential (through a 3 L). In 
the more realistic case of three-dimensional systems, this length-dependence translates into a volume-dependence^] 
of the harmonic force constants 

The volume dependence of the phonon frequencies induced by the volume-dependence of the force constants is 
traditionally modeled by the Griineisen parameter 

<91nzA (k) . 

™ = -^v 1 (24) 

which is defined for each branch j and each point k in the first Brillouin zone. But since we are interested in 
determining the free energy of a system, it is convenient to directly parametrize the volume dependence of the free 
energy itself. This dependence has two sources: the change in entropy due to the change in the phonon frequencies 
and the elastic energy change due to the expansion of the lattice: 

F (T, V) = E* (V) + F H (T, V) (25) 

where E* (V) is the energy of a motionless lattice constrained to remain at volume V , while Fjj (V) is the vibrational 
free energy of a harmonic system constrained to remain at volume V. The equilibrium volume V* (T) at temperature 
T is obtained by minimizing this quantity with respect to V. The resulting free energy at temperature T is then given 
by F(T,V* (T)).0 

Let us consider a particular case that illustrates the effect of temperature on the free energy, at the cost of a few 
reasonable assumptions. We assume that 

• the elastic energy of the motionless lattice is quadratic in volume; 

• the high temperature limit of the free energy can be used. 



Of course, in general, it could be a general strain dependence, if the symmetry of the crystal is sufficiently iow. 
16 Formafly, the free energy should be determined by a sum over every possible volume: — fcBTln \^2 V exp (— f3F (V))) • 
However, since the voiume is a macroscopic quantity, its distribution can be considered a delta function and the sum reduces 
to a single term: the free energy at the volume that minimizes the free energy. 
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As shown in Appendix [B[ in this approximation, the volume expansion AV as a function of temperature takes on 
a particularly simple form: 



AV _ 3k B Tj 
~W ~ B 



where 7 is an average Griineisen parameter: 



, 3N 

— Y 



V di 



3A ^ v m dV 

m—1 



(26) 



(27) 



The resulting temperature dependence of the free energy is given by 



F(T) _F (T, V ) (3k B Tj) 2 



N 



N 2B (V /N) ' 



(28) 



These expressions provide a simple way to account for thermal expansion. They also allow us to estimate the changes 
in vibrational entropy as a function of temperature that is due to thermal expansion: 



dS- 



vib 



dT B(V /N) 

In metallic alloys, this quantity is typically of the order of 10 _4 fc B /K 



(29) 



2. Simulation 



There arc two main simulation-based approaches to handling anharmonicity: Monte Carlo (MC) (Binder and 
Hccrmann, 1988) and Molecular Dynamics (MD) (Allen and Tildesley, 1987). While both approaches are able to 
model anharmonicity at any level of accuracy, they suffer from two limitations. First, they are computationally 
demanding and therefore have, to date, been limited to simple energy models. Second, they are unable to model 
quantum mechanical aspects of vibrations and are therefore limited to the high temperature limit There is an 
interes ting a nd useful comp lementarity between the quasi-harmonic model and simulation techniques flde Fontaine 



etaL, 1998; Morgan , 1998). Quantum effects typically become negligible in the temperature range where strong 



anharmonic effects, which cannot be modeled accurately within the quasiharmonic framework, become important. 

The use of sim ulat ion techniques to determine vibrational properties bypasses the coarse graining framework pre- 
sented in section [I C : Both configurational and vibrational excitations are treated on the same level. When a simple 
energy model provides a sufficient accuracy, one can calculate thermodynamic properties directly from MC simula- 
tions where both atomic displacements and changes in chemical species are allowed during the simulation (Rittncr 
et aL|, |1994|; |Silverman et al.|, |1995| ) . 



While a full determination of a phase diagram from simulations has, to our knowledge, not been attempted, both 
MD ( QRavelo et aT| , |1998| )) and MC ( Qdo Fontaine ct alj , |1998| ; |Morgan| , |199S| )) have been used to determine differences 
in vibrational free energy between two phases. Because neither MD nor MC are able to provide free energies directly, 
a special integration technique needs to be used. The idea is to express a thermodynamic quantity inaccessible to 
MC as an integral of a quantity that can be obtained through MC. A simple example is the change free energy F as 
a function of temperature at constant pressure, which can be derived from the Gibbs-Hclmholtz relation 



r l/T 2 

F{T 2 ) =F(Tx)- \ 

Jl/Ti 



Ed(l/T) 



(30) 



where E is the internal energy. Another example is the change in entropy as a function of temperature (at constant 
pressure) which can be expressed in terms of the heat capacity: 



S(T 2 ) = S(T 1 ) + £ ^dT. 



(31) 



7 Monte Carlo simulations that include quantum effects are possible for systems containing a small number of particles. 
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Often, the most computationally efficient path of integration between two states is not physically meaningful. For 
instance, one can gradually change the interatomic potentials during the course of the simulation, in order to model a 
change in the configuration of the alloy, without requiring atoms to jump between lattice sites. This task is achieved 
by defining an effective Hamiltonian 

H x = (1 - A) H a + XH 13 (32) 

that gradually switches from the Hamiltonian H a associated with phase a to the Hamiltonian Hp associated with phase 
(3 as the switching parameter A goes from to 1. This convenient path of integration permits the calculation of free 
energy differences between phases at a reasonable computational cost, with the help of the following thermodynamic 
relation: 

F? = F<*+ [ ({H ) x -(H a ) x )dX (33) 
Jo 

where (H a ) x is the average of the energy calculated using Hamiltonian H a (and similarly for (Hp)* . 



C. Energy models 



Force constants and anharmonic contributions arc ultimately always derived from an energy model. In this section, 
we discuss various energy models, from empirical potential models to first-principles techniques, and the error or bias 
they may introduce in the vibrational properties. 

Simple pairwise potentials of functionals (such as the Embedded Atom Method) are computationally efficient so 
that all vibrational properties can often be determined without any approximations beyond the ones associated with 
the specific energy model. For this reason, the use of simple energy models has proven to be an invaluable tool 



to understand trends in vibrational en t ropie s and to tes t a number of a pprox imations (Garbulsky and Ceder, 1996 



de Fontaine et al., 1998; Morgan et al, 1998; Garbulsky, 1996; Morgan, 1998) 



Several potential sources of error can arise when using pair potentials or pair functionals. The first one is that 
vibrational entropy is extremely sensitive to the precise nature of the relaxations that take place in an alloy and a 
simple energy model may not be able to accurately predict these relaxations. This problem is particularly apparent 
when considering the wide range of values found in the different calculations of the vibrat i onal entropy change upo n 



disordering of the Ni 3 Al compound (Ackland, 1994; Althoff et al., 1997; Ravelo et al., 1998; van de Walle et al., 1998) 



But, as shown in Table IV, most of the discrepancies can be explained from differences in the predicted volume change 
upon disordering. 

This is often aggra vated by the f act th at simple energy models are often not fitted to phonon properties. The 



problem was noted in ( Althoff et al , 1997) where the embedded atom potentials used were fitted to various structural 
energies and elastic constants (Voter and Chen, 1988). The acoustic modes were accurately extrapolated from the fit 



to the elastic constants, but the phonon frequencies associated with the optical modes were overestimated by about 
1O%.0 The question of the accuracy of simple energy models cl early merits furt h er at t ention. In this r espec t, the fit 



of simple energy models to the results of ab-initio calculations ( [Silverman et al. , 1995| ; Shaojun et al., 1998) offers a 
promising way to include vibrational effects. 

In oxides, electronic polarization has to be included in order to correctly model both the low frequency acoustic 
modes and the high frequency optical modes. Electronic polarization in oxides can be approximated with the so-called 



core and shell model (Dick and Ovcrhausser, 1958) 



While quantum mechanical methods are computationally more intensive, they generally provide more accurate force 
constants. The most obvious error introduced by the common Local Density Approximation (LDA) is its systematic 
underprediction of lattice constants which leads to an overestimation of elastic constants and phonon frequencies. This 
systematic error makes it difficult to compare the absolute values of calculated vibrational properties with experimental 
measurements. However, for the purpose of calculating phase diagrams, this bias may be less of a concern, because 
phase stability is determined by differences in free energies, and one would expect a large part of this systematic error 
to cancel out. 

A practical way to alleviate the LDA bias is to perform calculations at a negative pressure such that the cal c ulate d 
equilibrium volume agrees with the experimentally observed volume. As shown in (van de Walle and Ceder, 1999), 



Most of the bias in the vibrational entropy introduced by this problem should however cancel out when taking the difference 
in vibrational entropy between two phases where the same problem is present. 
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a very good estimate of the required negative pressure can be obtained by a concentration-weighted average of the 
pressure associated with the elemental solids. For the purpose of calculating elastic properties, this approach appears 
to outperform the most popular alternative to LDA, the Generalized Gradient Approximation (GGA).0 



D. Convergence issues 



The basic formalisms presented in Section II A and V A provide two natural ways to c ontrol the trade-off between 
accuracy and computational requirements. In the context of alloy theory (Section [I A ), the range of the effective 
clusters interactions included in the cluster expansion controls how accurately the configurational dependence of 
vibratio nal p roperties is modeled. In the context of the harmonic (or quasiharmonic) treatment of lattice vibrations 
(Section |V A| ) , the range of the force constants included in the Born- von Karman model controls the accuracy of the 
calculated vibrational properties for a given configuration. In principle, any desired accuracy can be reached, given 
sufficient computing power, by increasing the range of the interactions in both the cluster expansion and the Born- von 
Karman models. This section seeks to answer the important question of how far these two ranges of interactions need 
to be pushed in order to reach the accuracy required in a typical phase diagram calculation. 



1. Short-range force constant 



The evidence that spring models including only short-range force constants are able to correctly model vibrational 
quantities comes from various sources. 

First and second nearest neighbor spri ng model are routinely used to fit data obtained from neutron scattering 



measurement of phonon dispersion curves (Fultz et al., 1995; Nagel et al., 1995). In the theoretical literature, there have 



been direct studies of the convergence as a function of the range of interaction considered. All ab-initio studies find that 
short-range force cons tants (first or second neare s t neighbor) permit an accurate det ermination of thermodynamical 
quantities in met als ( van dc Wallc et al. , 1998 ; van de Wallc and G. Coder , 2000] ) and group IV semiconductors 
( parbulsky , 1996). It is important to note that this rapid convergence of most thermodynamic quantities occurs even 
when the pointwise convergence rate of the phonon DOS is slow. As noted before, this property arises from the fact 
that thermodynamic quantities are averages taken over all phonon modes and errors tend to average out. 

In ionic systems, the presence of long-range electrostatic interactions may require long-range force constants. How- 
ever, this electrostatic effect can easily be modeled using pair interactions at a moderate computational cost. Once the 
forces predicted from a simple electrostatic model have been subtracted, the residual forces should be parameterizable 
with a short-range spring model. 

Some of t he ab-initio studies of convergence hav e suggested additional simplifications to force constant tensors 
( |Garbulsk"y , 1996 ; van dc Wallc and G.Cedei , |2000 ): instead of attempting to compute all force constants in each 
tensor, is it possible to obtain reliable results by keeping only the largest terms. We now present a hierarchy of 
approximations that is a formalization of these findings. 

To obtain a more intuitive representation of a given force constant tensor <$> a p we express it in a basis such 

that the first cartesian axis is aligned along the line joining atom i and j. The second axis is then taken along 
the highest symmetry direction orthogonal to the first axis while the third axis is chosen so obtain a right handed 
orthogonal coordinate system. 

In the absence of symmetry, the most general force constant tensor has 9 independent elements. The first simplifi- 
cation, is to neglect the three body terms in the harmonic model of the energy (e.g. (x a (i) — x a (j)) (x@ (i) — (k)) 
with a 7^ 0). Physically, such terms arise from the deformation of the electronic cloud surrounding atom i that is 
caused by moving atom j and that affect the force acting on atom k. Clearly, for any force constant other than the 
nearest neighbor, this effect is negligibly small. Even for nearest neighbor tensors, it is the most natural contribution 
to neglect first. In can be readily shown that a solid consisting only of pairwise harmonic interaction, the tensor 
associated with a pair of atoms is symmetric: 



$a0 (i,j) = $0 a 

(This constraint is distinct from the conventional constraint: (i,j) — $/3 Q (J,i) ■ ) 



(34) 



1 Part of the success of the "negative pressure" LDA is due to the fact that it uses information regarding the true experimental 
volume instead of being fully ab-initio. But the knowledge of one pressure per element is a relatively small amount of information. 
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The elements of the force constant tensor can be ranked in decreasing order of expected magnitude based on three 
simple assumptions: 

1. Force constants associated with stretching a bond are larger than the ones associated with bending it. 

2. Terms relating orthogonal forces and displacements are smaller than those relating parallel forces and displace- 
ments. 

3. In the plane perpendicular to the bond, the anisotropy in the force constants is smaller than the magnitude of 
the force constants themselves. 



We then obtain 



with 



a d + e d — e 
= \ d + e b + c f 

d — e f b — c 



\a\ > \b\ > |c| > \d\ > \e\ > \f\. 



(35) 



(36) 



This hierarchy of force constants is important to keep in mind, given that the off-diagonal elements of the spring 
tensors are the most difficul t to obtain from superccll calcul ations, requiring much bigger supercells than diagonal 
elements. There is evidence (van de Walle and G.Ceder, 2000) that even keeping only the stretching (a) and isotropic 
bending (6) terms of the nearest neighbor spring tensor can provide vibrational entropies with an accuracy of about 
0.03 kB- If this observation turns out to be generally applicable, this offers a simple way to account for vibrational 
effects in phase diagram calculations. 



2. Short-range effective cluster interactions 



If a cluster expansion of the vibrational free energy only requires a small number of ECI to accurately model the 
configurational-dependence of the vibrational free energy, it then becomes practical to determine the values of these 
ECI from a small number of very accurate calculations of the vibrational free energy of a few structures. 

The issue of the speed of convergence of the cluster expansion is also related to the task of devising efficient ways 
to compute vibrational properties of disordered alloys: The faster the cluster expansion converges, the easier it is to 
model a disordered phase (see Appendix |C|). The calculations of the vibrational entro py change up o n disordering 



has proven to be a very effect i ve way to assess the importance of la ttice vibrations ( [Althoff ct al , 1997; Ravelo 



ct alj , |1998| ; [van de Walle ct alj , |1998| ; pvan de Walle and G.Ceder| , |200Cj) , since this quantity can be straightforwardly 



used to estimate the effect of lattice vibrations on transition temperatures with the help of Equation (^|) . 

The central question is thus whether the cluster expansion of the vibrational free energy converges quickly with 
respect to the number of ECI. This is a question distinct from the range of force constants needed to obtain accurate 
vibrational properties. The range of ECI needed to represent the conhgurational dependence of vibrational free energy 
may very well exceed the range of the force constants. Even in simple Born von Karman model s ystems, there is no 
direct correspondence between ECI and force constants, except in special case s (see Section VII A ). Once relaxations 



are introduced in the model, then all hope of a simple correspondence is lost (Morgan et al 



1998) 



In this context, the question of the existence of a rapidly converging cluster expansion of vibrational properties has 
to be answered through numerical experiments. Simple energy models offer the possibility to test, at a reasonable 
computational cost, the speed of convergence of a cluster e xpansion. Explicit calculation s of a well converged cluster 
expansion of vibrational entropy in a Lennard- Jones solid ( |Garbulsky and Cedci , |l996| ) have indicated that a small 
number of ECI (9) can provide a good accuracy (±0.03fcB). Other ben c hmarks of the speed of convergence , based 



on studies of disordered alloys ( ie Fontaine ct al. , 1998| ; Morgan et al. , 1998 ; Morgan , 1998 ; Morgan et al. , 2000 ) 



also indicate that concise and accurate cluster expansions are possible. Experiments tha t seek to link features of 



projected phonon DOS to the local chemical environment of the atoms ( Fultz et al. , 1998| ) suggest that short-range 



ECI should be able to successfully model vibrational entropy differences. 



One potenti al source of concern is the 

van de 



1998 



difficulty associated with accounting for t he siz e mismatch effect using a short-range ECI ( Morgan et al 
Walle and G.Cedei, 2000] : Morgan ct al., 200C). In the context of cluster expansions of the energy, relaxations of the 
atoms away from their ideal lattice site as a result of size mismatch are known to introduce both non negligible long 
range pair ECI and numerous multiplct ECI. A cluster expansion of the vibrational free energy is expected to exhibit 
the same problems. 
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All the full phonon DOS ab-initio calculation s of vibrational entropies in alloy systems performed so far have r elied on 
the ra pid convergence of the cluster expansion flGarbulskyj , |1996| ; |Tepesch et at] , |l996j ; rvan de Walle et all 1998r Dzolin: 



et afl, 1998c; van de Walle and G. Coder, 200C). While efforts to q u antify the error introduced by truncating the cluster 



expansion in ab-initio calculations have been made (Garbulsky , 1996; van de Walle et al 



1998 



van de Walle and 



G.C |ede^ , pOOOl ), the issue of the speed of convergence of the cluster expansion in the context of vib rational properties 



clea rly merits further study, especially in light of the importance of the size mismatch effect ( |van dc Walle and 



G.C|ed"er], |2000| ) 



VI. EXPERIMENTAL TECHNIQUES 



The experimental literature on the thermodynamics of lattice vibrations in alloys relies on mainly three techniques. 
In differential calorimetry measurements, the heat capacity of two samples in a different state of order is compared over 
a range of temperatures. If the upper limit of the range of temperatures is chosen to be sufficiently low, substitutional 
exchanges will not occur and the difference in heat capacity can be assumed to arise solely from vibrational effects. 
Integration of the difference in heat capacity (divided by temperature) then yields a direct measure of the vibrational 
entropy differences between the two samples of the range of temperature considered. This, of course, assumes that 
the lower temperature bound is sufficiently low, so that the vibrational entropy of both samples can be assumed to 
be zero at that temperature. It also assumes that the electronic contribution to the heat capacity is negligible. In 
practice, both assumptions are typically satisfied. The main problem with this method is that one is usually interested 
in vibrational entropy differences at the transition temperature of the alloy, which is usually above the upper limit of 
the temperature range used in the heat capacity measurements. The heat capacity therefore needs to be extrapolated 
to high temperature. This constitutes the main source of inaccuracies in this method. Examples of the use of this 



method can be found in ( 


Anthony et al., 


1993; 


Anthony et al., 


1994; 


Nagel et al., 


1997; 


Nagel et al., 


1995 


et al., 


L99I : 


Nagel, Fultz and Robertson 


, 1997). 



A second method is the measurement of phonon dispersion curve through inelastic neutron scattering. For ordered 
alloys that can be produced in large single crystals, this method is very powerful. Once the dispersion curves 
along special directions in reciprocal space are measured, they can be used to fit Born-von Karman spring models 
which, in turn, yield the normal frequencies for an y poi nt in the Brillouin zone. With the help of the standard 
statistical mechanics techniques described in Section VA, this information is sufficient to determine the vibrational 
entropy. Examples of ap plications of this method can be found in (Anthony et al., 1994; Nagel et al] 
' 1995) ; |Nagel et al.| , |1995| ) " 



et al, 1995; Fultz et al 



1997 



Fult; 



The applicability of this method is unfortunately limited by 
the availability of large single crystals. The case of disordered alloys presents an even more fundamental problem: 
Disordered alloys do not have well defined dispersion curves and there is no straightforward way to fit the spring 
constants of a spring model from the experimental data. This problem is usually addressed by using the virtual 
crystal approximation, in which different constituent atoms are replaced by one "average" type of atom (see Appendix 
^) . Unfortunately, this approximation has re peatedly been shown to have a very limi t ed accuracy for the purpose 



of measurin g vibrational entropy differences (Althoff et al 
|et al 



1997| ) 



1997; Ravelo et al, 1998; Shaojun et al., 1998; Nagel 



Nevertheless, single crystal phonon dispersion curve measurements for ordered alloys present a unique 
opportunity to perform a stringent test of the accuracy of theoretical models. 

A third method is the determination of the phonon density of states from incoherent neutron scattering measure- 
ments. In contrast to the preceding approach, this method can readi l y be applied to disordered systems and to 
compounds for which single crysta l s arc not availabl e ( Fultz et al. , 1995 ; Nagel et al. , 1996] ; Nagel, Fultz and Robert 



son, 1997; Bogdanoff et al., 1999) ; Robertson et al., 1999). The main limitation of this approach is that different 
atomic species have different neutron scattering cross-sections. The scattered intensity at each frequency measures a 
"density of states" , where each mode is weighted by the scattering intensity of the atoms participating in the mode 
in question. Thus, one needs some prior information about the vibrational modes in order to reconstruct the true 
phonon DOS from the experimental data. In the case of alloys, there is not a one-to-one correspondence between 
the measured data and the vibrational entro py. This problem can be alleviated by choosing alloy syste ms where the 



scattering intensity of each species is similar ( Nagel, Fultz and Robertson , 1997 ; Robertson et al. , 1999 ) 



Other techniques have been used to measure vibrational entropy differences. Some researchers have used the fact 
that vibrational entropy and thermal expans ion arc directly re l ated, to estimate vibrat i onal entropy differences from 

|1996| ; [Mukherjee et al.| 



accurate thermal expansion measurements (Mukherjee et al 



1998). The measurement of 



inelastic nuclear resonant scattering spectrum has also been used to relate changes in the phonon DOS to changes 
in the short-range order of a disordered alloy (Fultz et al., 1998). Finally, relatively noisy estimates of vibrational 



entropy differences can be obtained from X-ray Debye- Waller factors or from the measurement of mean square relative 
displacement (MSRD) of the atoms relative to their neighbors through extended electron energy-loss fine structure 
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(EXELFS) flAnthony et all |l993|) 



VII. MODELS OF LATTICE VIBRATIONS 

While the ability to control the level of approximation discussed in the previous section is extremely useful, there 
remains the problem that, very often, only considering the first few levels in this hierarchy of approximations already 
involves substantial computational requirements. For this reason, models of lattice vibrations that involve fewer 
parameters but more physical intuition may provide a practical mean of including vibrational effects in phase diagram 
calculations. In this section, we will present the advantages and weaknesses of each method, in light of the three 
fundamental mechanisms described in the Section |TV|. 



A. The "bond proportion" model 



There have been many at tempts (see, for instance, ( |Dyson| , |1953| ; |Wojtowicz and Kirkwoodj , |1960| ; |Bakkerj , |1982q ; 
Garbulsky and Ceder, 1994)) to find ways to express the relationship between the vibrational free energy and the 
dynamical matrix in a form that illustrates the intuition behind the "bond proportion" mechanism. In a variety of 
simple model systems, a convenient exact expression can be derived for the nearest neighbor ECI in the expansion of 
the vibrational free energy in the high temperature limit. For simple nearest- neighbor spring models with central forces 
in linear chains (Bakkei, 1982q; Matthew et aL, 1983 ;_ Garbulsky and Ceder, 1994), square lattices (Bakker, 1982c 



or simple cubic lattices (Waegmaekers and Bak 



ker 



1984 ) , the nearest neighbor ECI is given by 



d. rj,. ( k A A^BB 



(37) 



where kAAi ksB and \zab are, respectively, the spring constants associated with A — A, B — B and A — B bonds 
and d is the dimensionality of the system. It has been noted, on t he ba sis of numerical experiments, that the same 
expression performs well for other lattices (Garbulsky and Ceder, 1994). This success arises from the fact that, as 
shown in Appendix Equation ( |37j ) is the first order approximation to the true vibrational entropy change in a large 
class of systems which satisfies the following assumptions: 

• the high temperature limit of the vibrational entropy is appropriate; 

• the nearest-neighbor force constants can be written as <I>(i, j) = k ai<Tj cj>(i,j) where fc CTi(T denotes the (scalar) 
stiffness of the spring connecting sites i and j with occupations Cj and aj while the <p(i,j) are dimensionless 
spring constant tensors. The 4> ihj) are assumed equivalent under a symmetry operation of the space group of 
the parent lattice; 

• all force constants k n . a . are such that 



- 1 



< 1. 



(38) 



Equation ( |37| ) applies to simple harmonic models with nearest neighbor springs on the fee, bec or sc primitive 
lattices (and, approximately, on the hep lattice), as long as the above assumptions are satisfied. Both stretching 
and bending terms are allowed in the spring tensors, as long as their relative magnitude is independent of k aiaj (e.g. 
when the bending terms are always, say, 10% of the corresponding stretching term, regardless of the magnitude of 
the stretching term). If force constants, other than bending or stretching terms, are important, the bond proportion 
model ceases to be valid. This can be seen by the following argument. The "bond proportion" picture requires every 
bond of a certain type (for instance, A — A bonds) to have an identical spring tensor. However, the point symmetry 
of each bond can be different and similar chemical bonds in different environment face different symmetry-induced 
constraints on their spring tensors (Sluiter et al., 1999). The only way to reconcile these observations is use a spring 
tensor that is compatible with the highest possible symmetry, ensuring that it is also compatible with any other 
environment with a lower symmetry. With the highest possible symmetry, only two independent terms remain in the 
spring tensor: the stretching and bending terms. 

Equation ( |37| ) embodies the essential intuition behind the effect of the alloy's state of order on its vibrational free 
energy. When one replaces an A — A bond and a B — B bond by two A — B bonds, the vibrational free energy will 



23 



decrease only if the stiffness of A — B bonds, k^B, exceeds the geometrical average stiffness of the bonds between 
identical species \k~AAk~BB- This observation allows the determination of the expected effect of vibrations on the 
shape of the phase diagram by simple arguments. The link between the nearest neighbo r ECI of the expansion of the 
vibrational entropy can be summarized by the expression ( Garbulsky and Cedei , 1996 ): 



^config+vib 



T, 



con fig 



1 T aVi nn /k B T 



(39) 



where the "— " and "+" correspond to ordering and segregating systems, respectively, and where a is a dimensionless 
parameter that only depends on the lattice type and the ordering tendency of the system (for instance, for fee, a = 1.7 
in ordering systems and a = 9.8 in segregating systems, while for bec, a — 6.5 in both cases). 

It is straightforward to include vibrational effects in phase diagram calculations using the "bond proportion" model. 
All that is needed is an estimate of the stiffness of A — A, B — B and A — B bonds, which could come, for instance, 
from supercell calculations of the nearest neighbor force constants in a few simple structures or from the bulk moduli 
of the pure elements and one ordered compound. The nearest neighbor ECI then obtained can be simply added to 
the cluster expansion of the energy. 

While Equation ( |37j ) is useful to estimate the importance of the "bond proportion" mechanism in a given system, 
one can avoid some of the approximations involved in deriving Equation (37) at the expense of only a modest amount 
of additional effort. One can find the exact phonon DOS of the nearest neighbor Born-von Karman model for a 
variety of configurations of the alloy, which allows a more accurate cluster expansion of the vibrational energy to be 
derived. In this fashion, the condition specified in Equation J3lj| ) is no longer needed and the vibrational entropy can 
be calculated at any temperature. 

It is important to keep in mind that two important assumptions are made when invoking the "bond proportion" 
mechanism. First, vibrational entropies are solely determined by the nearest neighbor force constants. There is 
theoretical evidence that ne arest neighbor spring models can predict vibrational entr opy differences with an accuracy 
of ab out 0.02fce in metallic ( van de Walle et al. , 1998| ; van de Wallc and G. Coder . 2000 ) and semiconductor ( Garbulsky , 
1996) systems. Given that configurational entropy differences are typically of the order of 0.2fcs, this precision should 



be sufficient for practical phase diagram calculations. 

The second assumption is that each type of chemical bond is a ssumed to hav e an in trinsic stiffness that is indep endent 
of i ts environment. First -principles calculations on the Li-Al ( [gluiter et al. , 1996 ) and on the Pd-V system ( van dc 
Wal fe and G. Coder , 200C ) unfortunately indicate that the stiffness of a chemical bond does change substantially as a 
function of its environment. This problem is serious, as it considerably limits the applicability of the "bond proportion" 
model. These changes of the intrinsic stiffness of the bonds as a function of their environment are precisely the focus 
of the two other suggested sources of vibrational entropy changes. In summary, while the "bond proportion" model 
gives an elegant description of one of the mechanisms suggested to be at the origin of vibrational entropy differences, 
it completely ignores the two other mechanisms, namely, the volume and size mismatch effects. 



B. The Debye model 

Perhaps the most widespread approximation to the phonon DOS g(i>) is the Debye model (see, for instance, 
(Grimvall, 1986 ; Ashcroft and Mcrmin, 1976| )), where the phonon problem is solved in the acoustic limit. In this case, 
the phonon DOS is approximated by: 



9v 

= J IP- 



if v < vd 
if v > vr> 



(40) 



where vu = kB f L D and Qd is the Debye temperature, given by: 



Oz, - — 
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(41) 



where Cjj is the Debye sound velocity, defined by 



(42) 
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where the right-hand side is the directional average of a function of the three sound velocities C\ (Grimvall, 1986) 
The free energy of a Debye solid is given by: 

£- £ + |*.e D - W (D («?) - Sin (l - exp (-2£ 

« — — + 3fcsTln — ; — ) in the high temperature limit. 
N ' \ h J 



where the Debye function D (u) is given by 



D (u) = 3u 3 



T dx 



(43) 



Since the Debye sound velocity Co is a complicated function of all elastic constants of the material, an approximation 
to the Debye temperature that only involves the bulk modulus proves extremely useful. Such an approximation was 
derived by Moruzzi, Janak and Schwarz (MJS) (Moruzzi et al. , 1988) for cubic materials^: 



Q D = 0.617 



1/3 



k B V M 



1/2 



(44) 



where fl is the average atomic volu me B is the bulk modulus a nd M is the concentration weighted arithmetic mean 
of the atomic masses. As noted in (Garbulsky and Ceder, 1994), in the high temperature limit, the MJS model does 
not exhibit the required property that th e masses have no effect on th e vibrational free energy of formation, although 
using a geometric average of the masses (Garbulsky and Ceder, 1996| ) fixes this problem. 

The quasiharmonic approximation can be used, within Debye theory, to account for mild anharmonicity. In the so- 
called Debye- Griineisen approximation, the volume-dependence of the phonon DOS is modeled by a single Griineisen 
parameter and the effect of volume can be summarized by simply making the Debye temperature volume-dependent: 



Qd = Qd.o I ^ 



(45) 



where ©d,o is the Debye temperature at volume Vo and 7 is the griineisen parameter. 

Despite of its inaccurate description of the true phonon DOS at high frequencies, the Debye and Debye- Griineisen 
models are quite successful at modeling the changes in vibrational properties of a given compound as a function of 



temperature. For instance, the thermal properties of pure metals (Moruzzi et al., 1988) calculated in MJS approxi- 
mation are surprisingly accurate. The reason for this success is that most thermodynamic quantities (e.g. free energy, 
entropy, heat capacity, etc.) exhibit their most dramatic variations at low temperature, where the low frequency 
phonon modes that are correctly described by the Debye model have a dominant effect. In the high temperature 
regime, thermodynamic quantities are determined by the classical equipartition theorem, and any harmonic model 
gives the correct behavior. 

Debye-like models are expected to perform well in systems where the differences in vibrational free energy between 



( Ackland 




1994; 


Althofl 



vibrational e ffects in phase diagram cal c ulations and has res ulted in an improved agreement with experimental results 
( [Asta et al j |1993| ; [Sanchez et al.|, |199l|; |Colinet et al.| , |1994| ). 



However, as shown in (Garbulsky and Ceder, 1996] ), the Debye approximation and its successors can have significant 



shortcomings when used to calculate phase diagrams. A significant part of the vibrational free energy differences 
between different compounds arises from changes in the high frequency portion of the phonon DOS, which Debye-like 
models describe incorrectly. In some case s, the MJS approximation can even lead to an incorre ct prediction of the 



sign of the vibrational entropy difference (van de Walle and G.Cedei, 2000; Ozolins et al., 1998c) 



In summary, the Debye model and its derivatives capture the essential physics behind only one of the advocated 
mechanisms responsible for the configurational dependence of vibrational free energy: the volume effect. Approxima- 
tions based on the Debye model, however, fail to account for the possibility that the state of order also has a direct 
impact (i.e. not through the volume) on the shape of the phonon DOS (as predicted, for instance, by the "bond 
proportion" model). 



""although it has been used for materials with a lower symmetry (Asta et al, 1993) ; [Sanchez et al. , 1991) 



25 



C. The Einstein model 



The perf e ct co mplement to the Debye model is the Einstein model (see, for instance, ( Ashcroft and Mermin , 1976 
McQuarric, 1973| )), which focuses on the high frequency region of the phonon DOS, instead of its low frequency region 
The Einstein model assumes that a crystalline solid can be modelled by a collection of 3JV independent harmonic 
oscillators (3 per atom) sharing a common frequency. This frequency can, for instance, be determined by computing 
the natural frequency of oscillation of one atom when all others are frozen in place. This approach, known as the 
local harmonic model flLeSar et al. , 1989; Sutton, 1989|), has p roven especiall y useful to calculate vibrational entropies 
associated with defects^ Zhao et al. , 1993 ; [Najafabadi et al. , 1991 ; Sutton , 1989). The Einst ein model can also be 



combi ned with a Debye mode l to better fit experimental calorimctry data (Anthony et al., 1994) or thermal expansion 



data (Mukhcrjee et al., 1996) 



The local harmonic model is of little use whenever the system of interest exhibits translational symmetry, because 
the calculations required to determine the unknown parameters of an Einstein model from first-principles directly 
provide force constants. The latter could be used to obtain a more precise description of the DOS rather than the 
single-value DOS characterizing the Einstein model. 

The Einstein model is nevertheless extremely useful for conceptual purposes, as we will now illustrate. As shown 
in Appendix |d|, the vibrational free energy of a system is bounded by above and by below by the free energy of two 
Einstein-like systems: 



- 1 15=11* ((*-) 
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k B T 



(46) 



While the upper bound is obtained from the usual local harmonic model, where surrounding atoms do not relax, 
the lower bound is obtained when the surrounding atoms are allowed to relax freely. Another way to interpret these 
bounds is that at one extreme, each atom sees the others as having an infinite mass, while at the other extreme, 
each atom sees the other atoms as being massless. This result supports the view that vibrational free energy can be 
meaningfully considered as a measure of the average stiffness of each atom's local environment. 

A more rigorous way of definin g the contribut i on of an atom to the total vibrational free energy is the use of the 
projected DOS (see, for instance, ( Morgan et al. , 1998 )). This approach does not in any way simplify the calculation 



of vibrational properties, because the full phonon DOS is needed as an input, but it is a useful way to interpret the 
experimentally measured or calculated phonon DOS. To obtain the contribution of atom i to the DOS, the idea is to 
weight each normal mode by the magnitude of the vibration of atom i: 



9i W) 



(47) 



where ej is the eigenvector (normalized to unit length) associated with the mode of frequency Vj . Since extensive 
thermodynamic properties are linear in the DOS, atom-specific local thermodynamic properties can be readily defined 
from the projected DOS. Note that, by construction, all the projected DOS sum up to the true phonon DOS and 
thus, all the local extensive thermodynamic quantities sum up to the corresponding total quantity. 



D. The "bond stiffness vs. bond length" approach 



In ab-initio calculations, most of the computational burden comes from the calculation of the force constant tensors. 
It would thus be extremely helpful if the force constants determined in one structure could be used to predict force 
constants in another structure. From the failures of the "bond proportion" model, however, we know that forces 



constants obtained from one structure are not directly transferable to another structure (Sluiter et al 
Walle and G.Ceder, 



200C) 



199S ; van dc 



Nevertheless, a simple modification of the tran sferable force constant app r oach y ields substantial improvements in 
precision. First-principles calculation the Pd-V (van dc Walle and G.Ceder, 2000) system revealed that most of the 
variation in the stiffness of a given chemical bond across different structures can be explained by changes in bond 
length.^] This suggests that the transferable quantity to consider is a "bond stiffness vs. bond length" relationship. 



21 The results obtained in the Li-Al system (Sluiter et al, 199S| ) also suggest than length is a good predictor of stiffness, 
although these authors did not investigate this matter further. 
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As a first approximation, a linear relationship can be used 



a (I) = a a + ai (I — l ) 
6(0 =b + b 1 (l-l ) 

where a and b denote the stretching and isotropic bending terms, respectively and where ao and 60 describe the 
stiffness of the bond at its ideal length Zo while ai and bi are analogous to bond-specific Griineisen parameters. The 
other parameters of the spring tensor are unlikely to follow such simple relationships because they may be required 
to vanish according to the local symmetry of the bond, independently of its length (this is discussed in more details 
in Appendix ^) . 

This approximation was shown to be successful in the Pd-V system ( van de Wallc and G.Ced"er| , 2000| ). Figure \ 



illustrates the ability of this simple model to predict bond stiffness in diff e rent s tructures. A similar analysis performed 
with the data on the Ni-Al system from Reference (van de Walle et al. , 1998) is shown in Fig. |[ Table ^ compares 



the predictions obtained from the "bond stiffness vs. bond length" model with more accurate calculations. 

There are numerous advantages to this approach. From a conceptual point of view, this model presents a concise way 
to represent all three mechanism suggested to be the source of vibrational entropy differences. The "bond proportion" 
mechanism is the particular case obtained when little changes in bond length occur. The volume effect results from 
expanding all bonds by the same factor. The size mismatch effect (or the "stiff sphere" picture) is also modeled since 
the local change in stiffness resulting from locally compressed or expanded regions are explicitly taken into account. 
A straightforward way to represent the source of vibrational changes is to overlap the stiffness vs. length relationship 
and the changes in average length and stiffness in different states of order, as shown in Fig. ||. 

A second advantage of this method is computational efficiency. The unknown parameters of the model can be 
determined by a small number of supercell or linear response calculations. After that, the knowledge of the relaxed 
geometry of a structure is sufficient to determine the stiffness of all chemical bonds. Finding the vibrational entropy of 
the structure then just reduces to a computationally inexpensive Born-von Karman phonon problem. It is important 
to note that the knowledge of the relaxed geometries of a set of structures is a natural by-product of first-principles 
calculations of structural energies, which are needed to construct the cluster expansion of the energy in phase dia- 
gram calculations, whether vibrational effects are included or not. Since computational requirements do not grow 
rapidly with the number of structures considered, this opens the way for a much more accurate representation of the 
configurational-depcndcncc of the vibrational free energy. 

A third advantage of transferable bond stiffness vs. bond length relationships is that they contain all the information 
needed to account for thermal expansion as well, within the quasi-harmonic approximation. The slopes of the stiffness 
vs. length relationships for each chemical bonds explicitly defines the changes in phonon frequencies as volume changes. 
Since the bulk modulus of each structure is also a by-product of structural energy calculations, all the ingredients 
needed for a quasi-harmonic treatment of thermal expansion are available. 



VIII. CONCLUSION 



Lattice vibrations can have a significant impact on phase transition temperatures, short-range order, solubility 
limits, and the sequence in which phases appear as a function of temperature. The standard framework of alloy 
theory can be straightforwardly extended to account for lattice vibrations using the concept of coarse graining of the 
partition function. Once the degrees of freedom associated with lattice vibrations are integrated out, one is left with 
a standard Ising model, where the energy of each spin configuration is replaced by its vibrational free energy. The 
efficient evaluation of the vibrational free energy of each configuration is the main problem limiting the inclusion of 
lattice vibrations in phase diagram calculations. A number of investigations have sought to assess the importance of 
vibrational effects on phase stability, in order to ensure that the efforts involved in computing vibrational properties 
are justified. The conclusion of the most reliable of these studies is that vibrational entropy differences are typically 
on the order of 0.1 ks to 0.2 fee, which is comparable to the magnitude of configurational entropy differences (at most 
0.69 ks in binary alloys), thereby indicating that vibrations have a nonnegligible impact. 

The calculation of the vibrational free energy of a particular configuration of the alloy reduces to the well known 
phonon problem in crystals. While the standard harmonic treatment of this problem lacks the ability to model thermal 
expansion, which can have a significant impact on thermodynamic properties in alloys, this limitation is easily overcome 
with the help of the quasiharmonic model. An exact solution to the phonon problem for all possible configurations 
requires excessive computing power. However, the tradeoff between accuracy and computational requirements can be 
controlled in two ways, namely through the selection of the range of force constants in the Born-von Karman model, 
and through a choice of the number of ECI used to expand the configuration dependence of the vibrational free energy. 
While there is evidence that the range of force constants can be kept very small (first nearest neighbor springs), the 
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configurational dependence of the vibrational free energy is too complex to permit a drastic reduction in the number 
ofECI. 

Given the substantial computing power required to undertake lattice dynamics calculations, many attempts have 
been made to devise simpler models. For many years, the MJS approximation appeared to be a very promising way to 
include vibrational effects in phase diagram calculations, because it systematically improved the agreement between 
first-principles calculations and experimental measurements. This success may have been the result of two fortunate 
circumstances (i) First-principles phase diagram calculations typically overestimate transitions temperature and (ii) 
the MJS approximation nearly always yields a downward correction to the transition temperature. As the accuracy of 
phase diagram calculations improved th rough the use of longer-range cluster expansions (Tepesch et al., 1996; Ozolins] 
et aQ, [j|98 4 [van der Ven et alj , |l99q ), the systematic bias in the calculated transitions temperature substantially 
decreased. Simultaneously, more sophisticated models of lattice vibrations indicated that lattice vibrations do not 
always results in a reduction in the transition temperatures (Garbulsky , 1996; garbulsky and Ceder , 1996; van dc Walk 
and G. Coder, 2000). The net effect of these two trends is that, although the accuracy of first-principles calculations 
has increased over the years, obtaining improved agreement with experiment is now a much more stringent test. As a 
r esult, perfectly valid and accurate ca l culatio ns of vibrational effects sometimes reduce the agreement with experiments 



(Tepesch et al., 1996; Ozolins et al., 1998c). Hence, before one can unambiguously assess the importance of lattice 



vibrations through a full phase diagram calculation, all potential sources of error have to be carefully controlled, 
such as the precision of the energy model used and more importantly, the accuracy of the cluster expansion. To 
date, the most convincing evidence that taking lattice vibrations into account significantly improves agreement with 
experimental results comes from calculations of the l attice dynamics associated with a specific atomic configuration 
(e.g. a given compound or an isolated point defect) ( pzolins and Asta , 2001; Wolvcrton and Ozolins, [2001 ; Quong 



and |Lui| , 1997| ). In these settings, most sources of errors are under control and definite answers can be given 



Although the availability of more accurate computational tools has revealed that the trends in vibrational entropy 
differ ences between phases is far more complex that anticipated ten years ago ( ( Garbulsky and Cedei , [1996 



et aQ, [1997| ; [van de Walle and G.Cederj , |2000[ [Morgan et al| , [2000| ; [Bluitcr et all, , |1999| ; [Wolvcrton and Ozo insj , p00l| J)7a 
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simple picture of the mechanisms at work is now emerging. All the known so urces of vibrational entropy differences can 
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Morgan 



be con venie ntly summarized by the "bond stiffness vs. bond length" model (van de Walle and G.Ceder 
et aQ, |200C ). In this picture, each type of chemical bond is characterized by a length-dependent spring constant. 
Changes in vibrational entropy can originate from both changes in the proportion of each chemical bond and changes 
in their lengths as a result of local and global relaxations. This model not only provides an intuitive understanding 
of lattice vibrations in alloys, but also a practical way of including their effects in phase diagram calculations. This 
stiffness vs. length relationship of each type of chemical bond can be inferred from a small number of lattice dynamics 
calculations. The vibrational properties of any configuration can then be obtained at a very low computational 
cost from the knowledge of the equilibrium geometry of this configuration, an information that is already a natural 
by-product of any phase diagram calculation. 

Future investigations of the effect of lattice vibrations on phase stability should head towards three main directions. 

1. While reporting error bars is an important part of any experimentalists' work, theorists should devote signifi- 
cantly more effort to quantifying the uncertainties of their calculations. This would make it possible to clearly 
identify situations where the improved agreement with experimental results following the inclusion of vibrational 
effects is truly significant or merely the result of fortunate coincidences. It is admittedly difficult to quantify 
the errors introduced by the energy model (such as the LDA), but standard statistical techniques can clearly 
be used to quantify any error due to fitting the ab initio data with a simple model. 

2. Given the difficulty of extracting vibrational entropies from experimental data, theorists should undertake the 
computation of quantities that can be directly measured. For instance, a Born-von Karman model directly 
enables the simulation of incoherent neutron scattering data, while the inverse procedure is a highly non unique 
operation. The calculation of thermal expansion coefficients would also be a very sensitive test. 

3. There have so far been very few accurate phase diagram calculations that include the effect of lattice vibrations. 
The main limitation remains the determination of a cluster expansion that accurately models the configurational 
dependence of vibrational free energy. The "bond length vs. bond stiffness" model should prove to be an 
extremely useful tool in achieving this goal. Although this approximation has been very successful in all systems 
to which it has been applied, the confirmation of its validity in a wider range of systems is crucial. It would 
also be interesting to devise a hierachy of increasingly accurate approximations that would include the "bond 
length vs. bond stiffness" model as a particular case. 
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APPENDIX A: THE ABSENCE OF MASS EFFECTS IN THE HIGH-TEMPERATURE LIMIT 



This appendix shows that the vibrational entropy of formation is independent of the atomic masses in the high 
temperature limit, as several authors (Grimvall and Rosen, 1983; Garbulsky and Coder, 1996) have noted. In the 
high-temperature limit, the vibrational entropy is determined by the product of the frequencies of all normal modes 
of vibrations v m . which can be related to the eigenvalues A m of the 3iV x 37V dynamical matrix D of the system (up 
to a constant) 



x>("m) =in(n 

m \ m 




const 



2N 



In (det D) + const 



Using the properties of determinants, we can write: 

i In (det D)=^\n (det (m 1 ' 2 ^ 1 ' 

= ± l n (det ($) det (M)) = ± In (det ($) fj M*J 
= -^ln(det(*)) + A]T lnMi 



where M is the 3iV x 3iV diagonal matrix of all the N atomic masses of the system (each repeated three times) 
while $ is 3iV x 3N the matrix obtained by regrouping all the 3x3 force constant tensors $ (i, j) in a single matrix 
(analogously to Equation (|l4|)). Now consider the change in the value of ^2 m In {v m ) when an Na atoms of type A 
and Nb atoms of type B are combined to form an alloy. Let the subscripts AB, A and B respectively denote the 
properties of an A^ Na / n ^B^ Nb / n ^ alloy, a pure crystal of element A and a pure crystal of element B. 

A (y, In {v m ) ) =$> (vi B ) -x A ^ {vi) - x B J2 ln K) 

\ m / m m m 

= I In (det ($ AS )) + ^N A In M A + ^N B In M B 

-~ ln (det ($ A )) - ^N A In M A - \ In (dot - ^N B In M B 

= .- ln (dot ($ AB )) - - ln (dot - i ln (dot ($ B )) 

All the terms involving masses exactly cancel one another. 



To simplify the exposition and avoid the problem that the dynamical matrix has three zero eigenvalues associated with the 
possibility of a rigid translation of the system, we assume that some of the atoms of the system are attached to a fixed point 
of reference by a weak spring. In the thermodynamic limit, this assumption becomes inconsequential. 
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APPENDIX B: A SIMPLE MODEL OF ANHARMONICITY 



The material present ed in this appen dix combines standard results regarding the Griineisen framework that can be 
found, for instance, in (Grimvall, 1986). 

Two assumptions are made. First, the elastic energy of the motionless lattice is assumed quadratic in volume: 



E* (V) 



w^ 2 



(Bl) 



where B is the bulk modulus, Vq the equilibrium volume at OK (ignoring zero-point motion) and AV^ = V — Vq. 
Second, the high temperature limit of the vibrational free energy is used: 



F mb (T,V) = fcsT^ln 



k R T 



In this approximation, the volume-dependence of F vib takes on a particularly simple form: 

dF mb (T, V) _ 3Nk B Tj 



dV 



V 



where 



1 x - V ov l 



(B2) 
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(B4) 



is an average Griineisen parameter. In the high-temperature limit, an average Griineisen parameter can easily be 
defined, because the population of the phonon modes is no longer temperature-dependent, and any change in entropy 
can be unambiguously attributed to shifts in phonon frequencies. At lower temperatures, the changes in phonon 
population would need to be accounted for as well. 

If we assume that the volume-dependence of the vibrational free energy is linear in volume, we have: 



F (T, V) = E* (V) + F mb (T, V) 

= E* (V) + F mb (T, V Q ) + 

= ^-(AV) 2 +F mb (T, Vq) 

Minimizing this expression with respect to AV yields: 
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where 3fc ^ T7 is the coefficient of volumetric thermal expansion. The resulting temperature dependence of the free 
energy (for one given configuration a) is given by 



F(T) _ F{T,V Q ) 
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(B6) 



It is interesting to note that half of the vibrational free energy decrease due to thermal expansion is canceled by 
the energy increase of the motionless lattice. Hence, vibrational entropy differences originating from differences in 
thermal expansion between phases have, relative to other sources of vibrational entropy changes, half the effect on 
phase stability. 



APPENDIX C: MODELING THE DISORDERED STATE 



Although in phase diagram calculations, the use of the cluster expansion bypasses the problem of directly calculating 
the vibrational entropy of a disordered phase, there are cases where it is of interest to directly calculate the vibrational 
properties of the disordered state. For instance, in studies that seek to assess the importance of lattice vibrations 
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(Althoff et al. , 1997; Ravelo et al., 1998; van de Walle et al., 1998; Ozolins et al., 1998c), it is instructive to compute 



the vibrational entropy change upon disordering an alloy, since this quantity can be straightforwardly used to estimate 
the effect of lattice vibrations on transition temperatures with the help of Equation (Q) . Here are the most common 
methods used to model the disordered state. 

Perhaps the most obvious and brute force approach to modeling the disordered state is the use of a large supercel l 
wher e the occupatio n of ea c h site is chose n at random. This approach was chosen in all EAM calculations ( Ackland , 



1994 



it 



Althoff et al. , 1997; Ravelo et al., 1998; Morgan et al.| , 1998) as well as in other investigations (Shaojun 



1998). Unfortunately, it is generally not feasible in the case of ab initio calculations. 



et a 

The virtual crystal approximation (VCA) consists of replacing each atom in a disordered alloy by an "average" 
atom whose properties are determined by a concentration weighted average of the properties of the constituents. In 
the limit where the chemical species have nearly identical properties, this approximation is justified. This model 
has been commonly used to interpret neutron scattering measurements of phonon dispersion curves in the case of 
disordered alloys (Fultz et al. , 1995 ; Nagel et al. , 1995). It has also been used in a some theoretical investigations 
( |Cleri and Rosato , 1992; Pcrsson et al.| , 1999). However, the virtual crystal approximation has been repeatedly 
shown to be insufficiently accurate for the purp o se of c alculating diffe r ences in vibrationa l entro pies between distinct 



compounds (Althoff et al., 1997; Ravelo et al., 1998; shaojun et al., 1998; Nagel et al., 1997). Its weaknesses are 
numerous: It is unable to model "bond proportion" effects, volume effects and local relaxations. It also fails to give 

a mass-independent high temperature limit. 

Special Quasirandom Structures (SQS) (Zunger et al., 1990| ) combine the idea of cluster expansion with the use of 
supercells. SQS are the periodic structures that best approximate the disordered state in a unit cell of a given size. The 
quality of a SQS is described by the number of its correlations that match the ones of the true disordered state. There 
is thus a direct connection between cluster expansions and SQS: If there exists a truncated cluster expansion that is 
able to predict properties. of the disordered state there exists an SQS that provides an equally accurate description of 
the disordered state. 

SQS have been very success f ully u sed to obtain electronic and thermodynamic properties of disordered materials 



(see, for example, ( Hass et al. , 199C )).The accuracy of the SQS approach in the context of phonon calculations has 
been benchmarked using embedded atoms potential s which allow f or th e comparison with the "exact" 

" i~998l ) 



entropy of the disordered state with a large supercell ([Morgan et al. 



vibrational 

It has been found that, for the purpose of 
calculating vibrational properties, an SQS having only 8 atoms in its unit cell already provide a good approximation 
of the disordered state in the case of an fee alloy at concentration 3/4. While the performance of this small SQS is 
remarkable in a model system where local relaxations are disallowed, it tends to degrade somewhat when relaxations 
are allowed to take place. This effect can naturally be explained by the fact that relaxations are known to introduce 
important multibody terms in the cluster expansion, which translates into the requirement that the SQS must correctly 
reproduce the corresponding multibody correlations. 

The success of small SQS opened the way for the use of more accurate energy models to calculate vibrational 
properties of disordered alloys. SQS have been applied to the ab-initio calculation of vibrational entropy in disordered 



Ni 3 Al and Pd 3 V alloys (van de Walle et al., 1998; van de Walle and G.Ceder, 2000) 



APPENDIX D: THE EINSTEIN MODEL 



In the Einstein model of a solid, the free energy, in the high temperature limit, is given by 

k B T , ( h 



In det D 

2 \k B T 



2 \k B T 
k B T . ( h 



In — — det M det $ 

2 \k B T 

where D and $ are, respectively, the 3N x 3iV dynamical matrix and force constant matrix of the system while M is 
the matrix of the masses: 

My = SijMj. (Dl) 
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It can be shown (Van de Walk, 1995| ) that for any positive definite matrix (f> 

det$<IJ$ ri , (D2) 

i 

implying that 

F<M^n(-^nW) (D3) 



where the right-hand side expression is nothing but the free energy of the system in the Einstein approximation. A 
lower bound can be obtained by a similar technique, by using the inverse of the force constant matrix 

det<I>> (n(* -1 )«) • (° 4 ) 

The interpretation of the inverse <& is simple: It is the matrix that maps forces F exerted on the atoms to the resulting 
displacements u of the atoms. 

u = (D5) 

While $m is related to the oscillation frequency of a single atom when all other atoms are held in place, ..) 1 

is related to the oscillation frequency of an atom i when all surrounding atoms are allowed to relax so that the force 
exerted on them remains zero as atom i moves. Atom i has mass Mj while all other atoms are considered massless 
and relax instantaneously. Atoms located infinitely far away from atom i are held in place with an infinitesimal force. 

In conclusion, the free energy of a system is bounded by above and by below by the free energy of two Einstein-like 
systems: 



APPENDIX E: DERIVATION OF THE "BOND PROPORTION" MODEL 



This append i x generalizes the results found in ( Bakkei , 1982a ; Matthew et al. , 1985 ; Garbulsky and Ceder , 1994 
Bakker, 1982a; Waegmaekers and Bakkcr, 1984) in order to handle more general lattice types. We show that, in an 
important class of systems, the bond proportion model is in fact the first order approximation to the true change in 
vibrational entropy induced by a change in the proportion of the different types of chemical bonds. 

The alloy system is assumed to satisfy the following conditions: 

• the high temperature limit is appropriate; 

• the nearest- neighbor force constants can be written as = ef>(i,j) where denotes the (scalar) 
stiffness of the spring connecting sites i and j with occupations cr; and <7j while the 4>(i,j) are dimensionlcss 
spring constant tensors. The 4> ihj) are assumed equivalent under a symmetry operation of the space group of 
the parent lattice; 

• all force constants k„.„. are such that 



< 1. 



(El) 



Consider a <i-dimensional solid made of N atoms connected by springs of characterized by symmetrically equivalent 
tensors k<f> Without loss of generality, the masses of all atoms are set to unity since the formation entropies in 

the high temperature limit are independent of the atomic masses (see Appendix |Aj). In the high temperature limit, 
the vibrational free energy per atom is given by: 



k B T 
2N 



^ In A r , 



(E2) 
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where the sum is taken over the nonzero eigenvalues X m of the dynamical matrix D of the system. (The zero eigenvalues 
correspond the modes where the whole crystal moves rigidly. In the thermodynamic limit, these few missing degrees 
of freedom are inconsequential.) 

Because all springs in the system are equivalent to each other, matrix D can be written as 

D = kC (E3) 

where C is a matrix of dimensionless geometrical factors independent of fc but specific to the type of lattice. From 
this expression of D, it follows naturally that eigenvectors of D are independent of k and that its eigenvalues can be 
written as 

A m = kl m (E4) 

where the l m are geometric factors independent of k. 

Consider what happens to SVib when the stiffness of one of the springs is changed from k to k + Ak. Let AD denote 
the corresponding change in matrix D. To the first order, the resulting changes in the eigenvalues are given by: 

AA m = u' m ADu m (E5) 

where u m is the (dimensionless) eigenvector of D associated with eigenvalue A m . Since D is linear in the spring 
constants, we can write 

AD = Ak B (E6) 

where B is matrix of geometrical factors independent of k and Ak but specific to the type of lattice. While matrix B 
also depends on which spring is being modified, the matrices B corresponding to each spring are equivalent under a 
symmetry operation of the crystal's space group. The changes in the eigenvalues can then be expressed as: 

AAi = Aku' m Bu m 
= Ak g m 

where gi is a dimensionless number independent of k and Afc. 
Substituting these results into Equation ( |E2| ), we obtain: 

F vih = 5Z ln ( klm + Ak 9m) ■ (E7) 

m 

To the first order, we can express the vibrational entropy change as 

dF vib 



AF vib = 



dAk 



Ak 

Ak=0 



2N ^ kl r 



, 1 9m \ Ak 



k 

where G is a dimensionless geometrical factor depending only on the lattice type. 

In the limit of Afc -C fc, we can obtain the change in vibrational entropy due to a change in all the spring constants 
by simply summing the effect of the change Ak s in the stiffness of each spring s: 

AF vih = k B TGj2^r (E8) 

s 

To determine the value of G, we consider the following particular case for which the exact vibrational entropy change 
is known. Once the value of G is known, it can be used in any other case sharing a same lattice type. 
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In a solid bound by springs of stiffness k is given by, if the stiffness of all springs is increased by Afc, each eigenvalue 
A m becomes A m k+ ^ k and the vibrational entropy becomes: 

Vlb 2N ^ \ k 



2N ^ V k 

i v 

_ , k B T^ k + Ak 



fvib + fc s r^^ + o((Afc) 

iVd f ^Afc 
F vib + k B T- 



2NZN/2^ k 



= F vih 



ksTd Afc 



ZN ^ k 

s 

where Z is the number of nearest neighbors and ^ s denotes a sum over all nearest neighbor bonds. Since this result 



is exact to the first order, we can compare it to Equation (E£) and identify the unknown constant G to be -^jj. We 
thus obtain the following result: 

s 

We now turn to the problem of calculating the vibrational entropy of mixing in a binary alloy. We first define a 
normalized 3iV x 3iV dynamical matrix D as follows: 

(1,3) = — 7= = (E1U) 

where k ai<7i is the spring constant of an A — A bond if site i is occupied by a A atom similarly for a site occupied by 
a B atom. For the purpose of calculating free energy of formation, this normalized dynamical matrix gives the same 
result as the usual dynamical matrix because the factors in the denominator exactly cancel out, for the same reason 
masses cancel out (See Appendix |A|) . This transformation normalizes the spring constant associated with A — A bonds 
and B — B bonds to 1 while the spring constant associated with A — B bond becomes (kABj /V^AAkBs) where fcyts, 
kAA and ksB respectively denote the true spring constants of A — B, A — A and B — B bonds. The usefulness of this 



normalization is to extend the applicability of Equation (E9) to the case where kAA and ksB are very different 



Let us start with a phase separated mixture of A and B atoms. Let us think of this system as one where all atoms 
are identical but where the springs connecting them can be either one of three types A — A, B — B or A - B. Where 
the springs are placed defines which type of atom sits at each site. We now replace one A — A bond in the pure A 
phase by an A — B bond and one B — B bond in the pure B phase by an A — B bond. By Equation (p9|), the resulting 
change in vibrational entropy per atom is: 

AF ksTd / kAB ^ ^ab -A (Ell) 

ZN \y/kAAk~BB \fkAAk~BB J 



To satisfy the assumptions of the above derivation, we require that kAB I ykAAksB ^1- If we create a total number 
uab of A — B bonds, we perform the above operation hab/^ times and the vibrational entropy change is: 

A p l _ n_AB k B Td / k AB j , k AB A (m2) 

2 ZN \y/k A AkBB \fk.AAkBB J 



To the first order (when kAB I ykAAkBB < 1), this expression is equivalent to 

= UAB ksTd / k\ B \ 

N 2Z \kAAkBBj 1 ' 

The nearest neighbor ECI of the cluster expansion of the vibrational free energy is thus: 

v ^ = i kBT ^{^f-)- (E14) 
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APPENDIX F: INSTABILITY 



An extreme case of anharmonic ity occurs when t he en ergy s urface, in the neighborh o od of a configuration a, has 



no local minimum. As noted in ( Craievich et al. , 1997 ) and ( Craievich and Sanchez , 1997 ), this situation occurs 
sufficiently frequently to deserve a particular attention. A typical example of such a situation occurs when the fec- 
based LTo structure is unstable with respect to a deformation along the Bain path, which leads to a bec-based B2 
structure. While it is possible to construct a separate cluster expansion for the fee and bec phases, the fundamental 
question that arises is: What is the free energy of the Llo structure? Since it is unstable, the standard harmonic 
expression for the free energy can clearly not be used. 

One suggested solution to this problem, des cribe d in (Craievich and Sanchez, 1997), is to perform the coarse graining 
in a different order than presented in Section II C . The sum over configurations is performed first, and the vibrational 
properties of the configurational averaged alloy are then calculated. The main limitation of this approach is that it 
would be extremely difficult to compute the averaged vibrational properties by any other method than by the so-called 
virtual crystal approximation (see Section ^|). Another limitation is that it only addresses instabilities with respect 
to cell shape distortions, ignoring instabilities with respect to internal degrees of freedom (i.e. atomic positions). 

In this section, we pres ent another approach to solve the instability problem. We argue that the general formalism 
developed in Section |II C| can in fact be adapted to allow for instability. 

While the coarse graining technique is most naturally interpreted as integratin g out the "fa st" degrees of freedom 



(e.g. vibrations) before considering "slower" ones (e.g. configurational changes) (Cedei, 1993), the time scale of the 
various types of excitations is, in fact, irrelevant. The partition function is simply a sum over states which can be 
calculated in any order. As long as we can associate any vibrational state v of the system with a configuration cr, the 
coarse graining procedure remains valid. 

Under this point of view, it is clear that it does not matter whether there is even a local minimum of energy in 
the portion of phase space associated with configuration a. What is important, however, is that the neighborhood of 
configuration a in phase space is thoroughly sampled (i.e. that the constrained system is ergodic) over a macroscopic 
time scale. There is no need for ergodicity within the time scale of the configurational excitations. If the neighborhood 
of a given configuration a is not fully sampled before the alloy jumps to another configuration a' , it is still possible 
that the unsampled portion of phase space around a will be visited at a later time, when the system returns to the 
neighborhood of configuration a. The ergodicity requirement at the macroscopic time scale imposes the important 
but intuitively obvious constraint that the phase space neighborhood of configuration a cannot contain states that 
are associated to different phases of the system. 

This discussion shows that there is no fundamental limitation to the applicability of the standard coarse graining 
framework in the presence of instability. However, we still need to describe how the free energy of an unstable 
configuration could be determined in practice. The task is simplified by the fact that the free energy of an unstable 
stable does not need to be extremely accurately determined, because unstable states are relatively rarely visited, even 
at high temperatures. Nevertheless, it is important to assign a free energy to those unstable states, to ensure that 
the Ising model used to represent the alloy is well-defined. 

The free energy associated with one configuration can be obtained by integrating exp [— /3E(a, v)] with respect to 
v over the portion of phase space associated with a. In the classical limit, we can label the vibrational states v by 
the position each particle takes and the integration limits can be found by geometrical arguments. The quantum 
mechanical equivalent of this operation is complex,^] but unlikely to be needed in practice. The unstable states are 
essentially never visited at low temperatures, where a quantum mechanical treatment would be essential 0. 

Focusing on the classical limit, we consider an unstable configuration a. Let D be the dynamical matrix evaluated 
at the saddle point of the energy surface closest to the ideal undistorted configuration cr.f^] We consider that when 
the state v of the system is such that one atom moves away from its position at the saddle point by more than r, it 
should not longer be considered part of configuration a. For an instability with respect to internal degrees of freedom 



23 The quantum partition function can be written as the trace of the matrix exp (—f3H), where H is the (multibody) Hamilto- 
nian of the system. The trace can computed in any convenient basis and in particular one could use Dirac delta functions. In 
this fashion, it is possible to define a localized free energy by summing only over the delta functions located in the neighborhood 
of one configuration a. 

24 This observation is related to the fact that quasi-harmonic approximation, which allows a quantum-mechanical treatment, 
is accurate up to a temperature where the classical limit is reached. 

25 Is is possible that an unstable configuration a cannot be associated with a saddle point and the derivation would have to 
be modified. In particular the bounds of integration would have to be made asymmetric. 
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(atomic positions), a natural choice for r would be half the average nearest neighbor interatomic distance. For an 
instability with respect to unit cell deformation, r could be half the change in the average nearest neighbor distance 
induced by the displacive transformation. 

The boundedness of the portion of phase space associated with a can be accounted for by replacing the usual 
classical partition function associated with one normal mode of oscillation i by 

* ex P (-^ 2 ) <** y * (-^ a * s2 ) ds ( fi ) 

where Xi is the i-th eigenvalue of the dynamical matrix, h is Planck's constant and Li is a measure of the size of 
the phase space neighborhood of a along the direction associated with normal mode i. This size parameter can be 
expressed in terms of the parameter r just introduced. Let U{ (j) = e '.^ where ej is the i-th eigenvector of D and 




Mj is the mass of atom j. After normalizing u,- L so that J2j u i U) = ^> ^hc number of atom in the system, we can 
then write 

L t = r( max r \\u t (j) - Ul (F2) 

where the maximum is taken over all nearest neighbor pairs of atoms j,j'. This choice of integration bounds approx- 
imately defines a neighborhood of a such that no atom moves farther than r from its position at the saddle point 
(relative to its neighbors). In this approximation, the free energy of an unstable state is given by: 

(F3) 

where Vi is the frequency of normal mode i and where the error function for real or imaginary arguments is given by 

2?; C 1 22 

erf (u) = -== / e~ u s ds. (F4) 
V n Jo 

The suggested definition of the free energy of an unstable configuration has interesting properties. First, as the 
neighborhood size Li increases, the expression reduces to the usual harmonic expression. The effect of the correction 
is not limited to unstable modes: Modes that are so soft that it is likely that the motion of the atoms exceeds r are 
also affected. There may obviously be other definitions of Li. The above example simply gives an example of how it 
could be calculated. 

Going back to our initial example of the LIq — ► B2 instability, we can now outline how this problem could be 
handled within the traditional coarse graining scheme. Two separate clusters expansion need to be constructed, one 
for the bcc phases and one for the fee phases. But since we now know how to assign a free energy to the unstable Ll 
configuration, the fee cluster expansion can be successfully defined. The free energy attributed to the LIq configuration 
should be sufficiently high so that the free energy curve of the fee phase in the vicinity of 0.5 concentration will lie 
above the free energy curve of the bcc phase, as it should. The fact that both CVM or Monte Carlo calculations on 
the fee lattice would attribute a positive probability to Ll -like structures should not be regarded as a problem: This 
is precisely what will ensure that the calculated fee free energy curve lies above the bcc one. 

The discussion has so far been concerned with the expression of the partition function, which is the relevant quantity 
to consider when the phase diagram is calculated with the CVM or the low temperature expansion. Let us now consider 
the implications of this approach to Monte-Carlo simulations. Thermodynamic quantities derived from averages, such 
as the average energy, are obviously unaffected by the presence of unstable configurations. For quantities derived from 
fluctuations, such as the heat capacity, slight modifications are needed. In traditional Monte Carlo simulations, the 
heat capacity arising from vibrational degrees of freedom is consistently neglected, and any thermodynamic quantity 
obtained from Monte Carlo can be unambiguously interpreted as the configurational contribution. In the more general 
setting presented here, there is an overlap between vibrational and configurational fluctuations and the only way to 
obtain well defined thermodynamic quantities is to fully account for the vibrational fluctuations. Fortunately, there 
is a straightforward way to do so. The total variance of the energy (or any other quantity) can be exactly expressed 
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as a sum of the variance within each configuration a and the variance of the average energy of each configuration: 

<s 2 ) - (Ef = £ £ p„bL - (e E p ^ 

a -uGct \ er vEcr 




where P av is the probability of finding the system in state a,v while P a = X^e<x ^ >< ™ an< ^ = St>ecr p"^u»' The 
first term is the usual fluctuation obtained from Monte Carlo. The second term is a correction which takes the form 
of a simple configuration average of fluctuations within each configuration. The fluctuation of a system constrained 
to remain in the vicinity of configuration a is usually just as simple to determine as its average properties. In the 
case of energy, the fluctuations within each configuration are simply related the heat capacity of a harmonic solid. 

The main objective of this section was to show that there is no fundamental problem associated with unstable states 
in coarse graining formalism. While it is true that the free energy of an unstable configuration is not uniquely defined, 
once a particular way to coarse grain phase space is chosen, the free energy of all configurations can be defined in a 
consistent fashion. There are admittedly some practical issues to be resolved regarding the practical implementation 
of coarse graining in the presence of instabilities, but the approach suggested in this section indicates that these 
difficulties can be overcome. 
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(form): Vibrational entropy of formation from pure elements, 
rnd: Disordered solid solution. 
*: metastable compound 

f: calculated from the data presented in the paper. 

TABLE I. Calculated Vibrational Entropy Differences. See Table [n| for abbreviations. 



pot.: pair potentials rmonic approximation 

EAM: embedded atom method 

TB: tight-binding 

TB-LMTO: tight-binding LMTO 

LMTO: linear muffin-tin orbitals 

ASW: atomic spherical waves 

PP: pseudopotential calculations 

cal.: differential calorimetry measurements 
lxtal: single crystal phonon dispersion measurements 
INS: incoherent neutron scattering measurements 
anh.: anharmonicity included. 



QH: quasiharmonic approximation 

MC: Monte Carlo 

MD: molecular dynamics 

D: Debye model 

BvK: Born- von Karman spring model 

SC: supercell method 

LR: linear response 

VCA: virtual crystal approximation 

MJS: Moruzzi, Janak & Schwarz model 

CE: cluster expansion 

SQS: special quasirandom structures 



TABLE II. Abbreviations Used in Tables | and |m|. 
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TABLE III. Experimental Measurements of Vibrational Entropy Differences. See Table [lj| for abbreviations. 
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TABLE IV. Relation Between the Vibrational Entropy Change upon Disordering and the Volume Change upon Disordering 
in Various Theoretical Investigations of the M3AI Compounds 



Compound (Structure) 



"Stiffness vs. Length" Model 



Inn Spring Model 
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-5.57 
-5.54 



-4.39 
-4.47 
-4.54 
-5.55 
-5.57 



TABLE V. Comparison Between Vibrational Entropies Obtained from the 
from a First Nearest Neighbor Spring Model. 



'Bond Stiffness vs. Bond Length" Model and 
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e (a, v 



FIG. 1. The Coarse Graining Approach. The global energy surface of an alloy system, which gives the energy Ei of each state 
i, can be partitioned into a set of local energy surfaces, each of which is associated with a distinct configuration Ok of an Ising 
model. Within each configuration , the local energy surface gives the energy as a function of the atomic displacements v (or 
any other nonconfigurational degrees of freedom). The thermodynamics of the system thus reduces to the one of an traditional 
Ising model, where the "energy" of each configuration now becomes the free energy associated with the corresponding local 
energy surface. If the local energy surfaces can be considered quadratic, their associated free energy can be determined by 
solving a Born-von Karman phonon problem. 




FIG. 2. The "Bond Proportion" Mechanism. Ordered alloys are characterized by the fact that a large proportion of the 
nearest-neighbor bonds join unlike atoms. As these types of bonds are presumably stiffer, they are responsible for an increased 
density of high-frequency optical phonon modes. Disordering reduces the proportion of stiff nearest-neighbor bonds, and the 
density of high frequency optical modes decreases correspondingly, resulting in an increase in vibrational entropy. 
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Ordere 




Disordered 



FIG. 3. The Volume Mechanism. Disordering is typically associated with an increase in volume. Since chemical bonds tend 
to soften as they lengthen, an overall increase in volume should be associated with a corresponding decrease in the frequency 
of all phonon modes, resulting in an increase in vibrational entropy. 
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FIG. 4. Stretching (s) and Bending (b) Terms of the Nearest-Neighbor Spring Tensor as a Function of Bond Length. Each 
point corresponds to one type of bond in one of a set of fee structures (LI2, DO22, SQS-8, fee Pd and fee V, each taken at two 
different volumes). 
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2.3 2.4 2.5 2.6 2.7 2.8 2.9 



FIG. 5. Bond Stiffness as a Function of Bond Length in the Ni-Al System. Only stretching terms are shown. Solid lines are 
the result of a fit to a second order polynomial. Each point corresponds to one type of bond in one of a set of fee structures 
(LI2, DO22, SQS-8, fee Al and fee Ni, each taken at two different volumes) 
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FIG. 6. Shift in Average Bond Stiffness (Along the Stretching Direction) and Bond Length upon Disordering. The fitted 
line of Fig. H is shown for reference. 
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